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efficiency  allowed  the  analysis  of  much  longer  filters. 
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Abstract 


A  variety  of  approaches  to  the  digital  filtering  of  voice 
signals  corrupted  by  background  engine  noise  are  presented. 

These  approaches  include  the  standard  Least  Mean  Squares  (LMS) 
adaptive  noise  cancellation  algorithm,  an  optimum  fixed  weight 
filter,  a  special  type  of  notch  filter,  and  frequency  domain 
adaptive  noise  cancellation.  The  filters  have  been  Implemented 
both  off-line  using  FORTRAN  programs  on  an  LSI-1 1/2  microcomputer 
and  in  real  time  using  the  Texas  Instruments  TMS  320 
microprocessor  and  in  PDP-11  assembly  language  using  an  LSI-1 1/2. 

Frequency  domain  adaptive  filtering  was  seen  to  be  superior 
to  the  LMS  time  domain  algorithm  because  it's  much  greater 
computational  efficiency  allowed  the  analysis  of  much  longer 
filters. 

A  digital  filter  that  exploits  the  periodicity  of  the  engine 
noise  by  having  a  notch  at  each  harmonic  was  seen  to  be  more 
effective  than  any  of  the  adaptive  filters.  The  digital  filter 
equations  are  derived  starting  with  an  adaptive,  finite  impulse 
response  filter  with  a  particular  periodic  reference  input.  This 
adaptive  filter  is  shown  to  be  in  reality  an  infinite  impulse 
response,  time  invariant  filter. 
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I.  Introduction 

Background  noise  at  the  source  has  long  been  a  problem 
affecting  the  intelligibility  of  Coast  Guard  voice 
communications.  Recent  studies  of  other  researchers  [7-10] 
indicate  that  additive  noise  will  become  a  much  more  severe 
problem  as  the  Coast  Guard  converts  from  analog  systems  to  voice 
privacy  systems  using  speech  compression  and  digital 
transmission.  In  particular,  Gold  and  Tierney  at  MIT's  Lincoln 
Laboratory  [7]  have  found  that  computer  simulated  F-15  cockpit 
noise  added  to  clean  speech  reduces  diagnostic  rhyme  test  (DRT) 
[11]  intelligibility  scores  to  92.6%  without  speech  compression 
but  to  75.2%  when  processed  through  a  2400  bit/second  LPC-10  [12] 
speech  compression  algorithm. 

8ur.  efforts  for  the  past  eighteen  months  have  focused  on 
determining  if  digital  filtering  in  general  and  adaptive  noise 
cancellation  in  particular  can  improve  the  signal  to  noise  ratio 
at  the  input  to  the  system  and  thus  the  intelligibility  of  the 
received  signal.  A  wide  variety  of  approaches  have  been 
attempted  with  varying  degrees  of  success.  A  summary  of  these 
efforts  is  contained  in  the  main  body  of  this  report.  The 
specific  details  including  program  listings  are  in  Appendixes  A- 
F. 


II.  Experimental  Methods 

The  general  technique  employed  has  been  to  make  stereo  tape 
recordings  of  background  engine  noise  both  with  and  without  a 
voice  signal  on  one  of  the  channels,  and  then  to  process  and 
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in  the  laboratory.  Recordings ' were  made 
my  Luder  yawl  with  the  engine  (a  four 
el)  running,  in  the  Academy  Power 
ith  a  single  cylinder,  four  cycle, 

7  running  and  most  recently  in  the  Power 
ith  a  four  cylinder  Datsun  710  running. 

:D  reel  to  reel  tape  deck  donated  by  CDR 
for  the  recordings.  All  recordings  were 
;nd.  General  Radio  Model  1952  Universal 
j  pass  anti-aliasing  and  reconstruction 

■ssing  the  recorded  waveforms  were  of  two 
gments  (less  than  one  second)  were 
;mory  .using  a  Digital  Equipment 
.rocomputer  and  a  Data  Translation  2785 
’liters  were  then  implemented  and  the 
le  using  FORTRAN  programs, 
is  to  implement  the  filter  in  real  time 
m  an  oscilloscope  or  to  listen  to  the 
i  one  exception  this  was  done  using  a 
JO  Evaluation  Module  (E7M).  The  notch 
indix  D  was  implemented  in  real  time  on 
the  TMS  320  EVM.  In  our  proposal  [133 
real  time  implementation  were  proposed. 

:W  TDC-1010  hardware  multiplier/ 

;ed  to  the  LSI-11/2  computer  but  even  with 
■  slow  to  implement  reasonable  size 
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adaptive  filters  in  real  time.  It  was  hoped  to  play  the  tape  - 
deck  into  the  filter  at  a  much  slower  speed  than  it  was  recorded, 
tape  record  the  filter  output  and  then  play  the  output  back  at 
the  original  speed,  but  we  were  unable  to  obtain  a  tape  deck  with 
this  capability.  Mo  attempt  was  made  to  use  the  Western  Digital 
WP  3150  Programmable  Digital  Filter.  We  were  unable  to  obtain  a 
final  data  sheet  and  have  concluded  the  announcement  was  a  trial 
balloon  and  that  the  circuit  never  went  into  mass  production.  An 
overwhelming  conclusion  we  have  reached  in  the  context  of  our 
work  is  that  the  TMS  320  is  by  far  the  best  approach  to  audio 
frequency  digital  signal  processing.  The  second  generation  of 
the  circuit  [14]  was  announced  in  February  1985  at  the  IEEE 
International  Solid  State  Circuits  Conference  and  will  make  even 

more  sophisticated  professing  possible; 

. 

III.  Summary  of  Approaches  and  Results 

a.  Analysis  of  LMS  Adaptive  Filters  in  FORTRAN 

Our  first  approach  was  an  off-line  FORTRAN  implementation  of 
the  standard  LMS  adaptive  noise  cancellation  algorithm  [1], 
Several  such  programs  were  written,  one  of  which  is  described  in 
Appendix  A.  When  the  inputs  were  highly  correlated  signals  from 
laboratory  signal  generators  the  filter  was  seen  to  be  very 
effective.  When  the  inputs  were  two  microphone  recordings  of 
engine  noise  there  was  virtually  no  improvement  in  signal  to 
noise  ratio. '  It  is  felt  this  was  due  to  a  combination  of  two 
problems.  The  filter  length  had  to  be  kept  short  to  insure  the 
weights  would  converge  within  the  sample  size  allowed  by  the  56 
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Kbyte  random  access  memory  of  the  LSI-11/2.  This  resulted  in 

filter  lengths  much  shorter  than  the  fundamental  period  of  the 

engine  noise  and  very  poor  filtering.  When  the  filter  length  was 

increased  the  weights  would  not  converge  within  the  allowable  -  ‘s 

sample  size.  The  recent  acquisition  of  an  LSI-11/23 

microcomputer  with  248  Kbytes  of  random  access  memory  from  the 

USCG  Electronics  Engineering  Laboratory  will  permit  the  analysis 

of  longer  filters  in  the  future. 

b.  Analysis  of  Optimum  Fixed  Weight  Filters 
It  was  next  attempted  to  eliminate  the  slow  weight 
convergence  problem  by  calculating  the  optimum  (in  a  linear  least 
squares  sense)  fixed  weight  filter  for  the  particular  data  set 
and  then  calculating  the  output  using  these  fixed  weights.  A 
.  description  and  listing  of  the  FORTRAN  program  to  accomplish  this’  , ,  , 

is  contained  in  Appendix  B.  Improvements  in  signal  to  noise 
ratio  of  6-8db  were  possible  when  the  sampling  rate  was  reduced 
to  1-2  kHz  and  the  anti-aliasing  filters  adjusted  accordingly,. 

This  effectively  increases  the  filter  length  in  time  to 
approximately  the  period  of  the  noise.  Although  these  low 
sampling  rates  are  unacceptable  for  voice  communications  the 

■« 

intent  was  to  determine  if  longer  filters  would  be  effective  at 

the  necessary  higher  sampling  frequencies.  It  is  also  believed 

there  were  numerical  problems  in  the  Gauss  elimination  subroutine 

when  attempting  to  analyze  longer  filters.  The  subroutine  should  ;.■! 

have  used  double  precision  arithmetic  [15]  but  single  precision 

was  used  due  to  limited  memory  and  because  double  precision 

floating  point  instructions  are  not  part  of  the  assembly  language 

•  -.-‘i 

"I 
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instruction  set  on  the  LSI-11/2  [18]  and  would  have  been  slow  to 
execute.  The  LSI-11/23  with  it's  larger  memory  and  double 
precision  instructions  in  assembly  language  will  allow  analysis 
of  longer  filters. 

c.  Real  Tiae  Adaptive  Filtering  Using  the  TMS  320 

At  this  point  (November  1983)  a  TMS  320  Evaluation  Module  was 
procured  and  efforts  shifted  to  real  time  implementation.  The 
first  version  was  done  by  then  Cadet  1/c  now  ENS  M.  D.  SAKAHARA 
as  his  Systems  Design  and  Synthesis  project.  His  report  was 
included  in  the  July  1984  interim  report  [16].  His  filter  had  40 
weights  and  would  operate  at  sampling  frequencies  of  up  to  8.5 
kHz. 

A  second  version  (Appendix  C)  used  BASIC  pn  the  Dartmouth 
Time  Sharing  System  to  write  much  more  efficient  TMS  320  code  and 
resulted  in  a  68  weight  filter  that  would  operate  at  sampling 
rates  of  up  to  10.7  kHz.  The  essential  difference  is  that  all 
the  looping  is  done  in  BASIC  at  the  assembly  language  generation 
level  resulting  in  straight  line  code  that  executes  at  the 
maximum  benchmark  speed  of  1.2  microseconds/weight  [17].  The 
first  version  used  conditional  branching  and  executed  at  2.6 
microseconds/ weight. 

Extensive  tests  were  done  on  recorded  engine  noise  and  the 
results  are  contained  in  Appendix  C.  The  conclusion  reached  was 
that  68  weights  are  not  enough  for  an  8  kHz  sampling  frequency. 
The  second  generation  TMS  32020  [14]  has  544  words  of  on  chip 
data  memory  compared  to  144  for  the  TMS  32010,  efficient  access 
to  64K  words  of  external  memory,  and  has  combined  some  operations 


that  were  two  instructions  into  a  single  200  ns  instruction. 
Adaptive  filters  of  at  least  120  weights  at  8  kHz  should  be 
possible. 

d.  Notch  Filters 

The  highly  periodic  nature  of  the  engine  noise  implies  that 
any  reference  input  that  contains  all  the  harmonics  of  the  noise 
in  the  primary  input  should  suffice.  A  reference  input  of  unit 
impulses  periodic  at  the  same  frequency  as  the  engine  noise 
contains  these  harmonics  and  results  in  very  efficient  algorithm 
allowing  filters  of  virtually  arbitrary  length.  Synchronization 
is  obtained  by  triggering  the  A/D  converter  with  an  optical  shaft 
encoder  mounted  on  the  engine. 

Further  analysis  of  the  algorithm  revealed  that  what  started 
as  an  adaptive,  finite  impulse  filter  was  in  fact  a  fixed  ‘weight, 
infinite  impulse  response  filter.  Preliminary  results  of  using 
the  filter  were  contained  in  the  interim  report  [16]  and  a 
detailed  analysis  is  contained  in  Appendix  D. 

This  filter  appears  to  be  the  most  promising  at  this  stage. 

It  is  intended  to  further  evaluate  it's  performance  by  making  a 
series  of  comparative  diagnostic  rhyme  tests  (DRT's)  [11]. 
Comparisons  of  intelligibility  will  be  made  among  filtered  and 
unfiltered  speech  using  both  an  omnidirectional  microphone  and  a 
Shure  Model  562  noise  canceling  microphone  [19=20].  Due  to  the 
failure  of  the  motor  generator  set  and  therefore  the  lack  of  DC 
power  it  has  not  been  possible  to  start  the  engines  in  the  Power 
Engineering  Laboratory  since  September  1984.  We  have  designed 
and  constructed  a  high  power  DC  supply  and  hope  to  have  engines 
running  in  the  near  future. 


The  LTiS  filter  was  inolenented  through  simulation  on  the 
LSI-11/?  computer.  Since  the  LSI-11  could  not _  execute  the 
adaptive  algorithm  fast  enough  for  real  time  applications,  data 
sets  '/ere  obtained  and  the  adaptive  algorithm  apnlied  to  the  data 
sets . 

A  Tortran  IV  program  './as  written  to  implement  the  adaptive 
filter  and  used  assembly  language  subroutines  to  communicate  with 
the  real  world  (  see  flow  chart  ) .  The  input  data  set  was  placed 
into  two  inout  arrays,  one  for  the  samples  with  noise  (  I.T1  )  and 
the  other  for  samnles  containing  noise  plus  the  desired  signal  ( 
I'.:?  ) . 

The  program  first  implements  a  fixed  weight  filter  on  the 
reference  inout  and  subtracts  the  resuLt  from  the  primary  input. 
The  output  of  the  adaptive  filter  was  pLotted  and  displayed  on 
the  video  terminal.  This  result  was  applied  to  the  adaptive 
algorithm  to  adjust  the  weight  vector.  Once  the  adaptation  cycle 
is  connle ted ,  the  computer  then  gets  the  next  input  and  starts 
the  process  all  over. 

'./hen  all  of  the  primary  samples  are  u3ed,  the  program  plots 
the  innulse  response  of  the  filter  and  calculates  the  frequency 
soectrum  if  the  weights.  noth  phase  and  magnitude  were  plotted. 
Since  the  adaptive  cycle  was  applied  after  the  fixed  filter,  the 
plotted  weight  vector  has  gone  through  an  extra  adaptation  cycle. 

In  the  strictest  case,  the  performance  of  the  plotted 
weights  is  unknown.  However,  this  dose  not  pose  a  “significant* 
problem  since  the  weights  change  very  little  once  they  have 
converged.  This  is  evident  in  the  plots  obtained  from  tne 
program  (  see  Appendix  1  ) .  Once  the  weights  converged  on  the 
L/-1S  solution,  they  change  very  little. 

The  entry  of  pertinent  constants  is  executed  at  the 
beginning  of  each  run.  On  the  first  run,  the  input  channels  ( 
17C1  and  'TO?  )  were  selected  and  the  number  of  samples  (  1TUV.S  )  , 
delays  (  1TUMD  )  and  weights  (  JJUMtf  )  for  the  filter  entered.  The 
adaptive  coefficient  (  ADAPT  )  was  entered  and  the  weight  ( 
'rsiGHT  )  and  input  arrays  were  reset  as  necessary.  The  data  entry 
segment  of  the  program  was  placed  at  the  and  of  the  program  to 
halo  speed  execution,  however  it  did  not  have  a  significant 
effect . 


no lamentation 


The  least  mean  square  (  L'**3  )  algorithm  of  the  adaptive 

filter  was  used  since  it  was  one  of  the  more  simole  forms.  It 
has  the  form  of 

•/C<+1)  »  ’■!(<)  +  u  Y  U(k) 

where  J(k)  is  the  weight  vector,  Y  is  the  output  of  the  adaptive 
filter,  U(k)  is  the  reference  input  array  and*  u  is  the  adaptive 
coefficient. 

The  reference  signal  was  filtered  and  then  subtracted  from 
the  primary  signal.  The  result  is  treated  as  an  error  signal  for 
the  filter.  For  a  positive  output,  the  adaptive  algorithm  will 
tend  to  increase  the  magnitude  of  the  weight  vector.  The  output 
of  the  filter  will  increase  causing  the  output  of  the  adaptive 
filter  to  decrease  or  tend  negatively.  'Then  the  output  goes 
negative,  the  algorithm  will  again  tend  to  alter  the  weights  to 
decrease  the  output.  This  tendency  to  minimize  the  output  is 
where  the  L»'tS  algorithm  gats  its  name. 

It  would  seem  that  if  this  filter  were  working  in  the  ideal 
case,  the  weights  would  converge  to  completely  cancel  the  output. 
The  adaptive  coefficient  is  what  prevents  this  from  happening. 

•The  magnitude  of  the  adaptive  coefficient  must  be  small..  Since 
the  adaptive  coefficient  has  a.  direct  affect  on  the  rate  of 
convergence  of  the  weights,  a  small  value  for  u  would  reduce  the 
tsndency  to  completely  cancel  the  output  and  yet  suppress  most  of 
the  interference. 

The  cancellation  of  the  noise  in  the  primary  input  is 
accomplished  by  reshaping  the  reference  input,  through  filtering, 
to  reflect  the  noise  in  the  primary  input.  If  the  filter  Is  to 
adjust  its  characteristic  to  accomplish  this  reshaping,  it  would 
be  desirable  for  the  characteristic  to  change  slowly.  This  avoids 
oscillating  around  the  desired  response  as  the  case  would  be  if 
the  weights  were  allowed  to  converge  too  quickly.  \  large 
adaptive  coefficient  may  reduce  the  effectiveness  of  the  overall 
filter  to  the  >  point  of  uselessness  (  greater  than  2.3*13  *9  in 
this  application  )  or  cause  the  filter  to  become  unstable  and 
blow-uo  (  u  greater  than  2.3  ).  <\  coefficient  that  is  too  small 

will  impede  the  weights  from  converging  or  even  completely 
preventing  them  from  reducing  any  noise  in  the  primary  input.  A 
compromise  must  be  reached  between  rate  of  convergence  and 
accuracy  of  the  convergence  in  the  proper  selection  of  u.  The 
adaptive  coefficient  must  be  small  enough  to  allow  the  filter  to 
accurately  reshape  the  reference  input  and  allow  the  weights  to 
converge  in  a  reasonable  time. 
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Introduction 


A  connor.  method  of  estimating  a  voice  3ignal  distorted  by 
additive  noise  is  to  pass  it  through  a  filter  that  tends  to 
reduce  the  noise  level  in  the  signal.  This  filter  can  either  be 
fixed  or  adaptive.  The  fixed  filter  requires  a  prior  knowledge 
of  the  signal  end  noise  in  order  for  it  to  be  effective. 

The  adaptive  filter  does  not  require  any  or  little  knowledge 
of  the  signal  or  noise.  It  has  the  ability  to  alter  its  impulse 
response  automatically  to  suppress  the  noise  fron  the  signal. 
The  adaptive  algorithm  used  in  this  application  has  two  inputs*  a 
primary  input  that  contains  the  signal  (  s.  )  and  noise  (  n,  ) 
and  a  reference  input  that  contains  noise  (  n*  )  that  is  somehow 
correlated  with  the  noise  (  n»  )  in  the  primary  signal  The 
reference  input  should  contain  none  or  very  little  of  the  signal 
since  the  filter  may  converge  to  cancel  to  signal  as  well  as  the 
noise  in  the  primary  input. 

\ 


l3iOCK  Ol  ap  PIUTStt. 


The  adaptive  filter  can  be  used  to  eliminate  the  5-T  llz  hum 
in  EKG's.  The  lldvac  outlet  can  be  used  as  the  reference  input 
and  the  normal  chest  leads  as  the  primary  incut.  In  a  similar 
application,  the  mother's  heartbeat  can  be  filtered  away  from  an 
infant' 3  EKG.  The  reference  leads  are  placed  near  the  mother's 
heart  and  the  primary  leads  are  placed  on  the  mother’s  abdomen. 
Interference  entering  the  side  lobes  of  a  receiving  antenna  can 
be  suppressed  by  using  a  circular  array  of  primary  antennas  and  a 
reference  antenna  chosen  in  the  area  of  the  incoming 
interference.  In  voice  communications,  where  background  noise  is 
a  problem,  two  microphones  may  be  used.  The  operator  speaks  into 
the  primary  microphone  while  the  reference  microphone  is  placs  in 
the  relative  vicinity  of  the  operator  picking  up  the  background 
noise . 
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Abstract 


An  LMS  adaptive  filter  was  implemented  on  the  L3I-11/2  computer, 
"he  computer  could  not  process  fast  enough  for  a  real  time  application 
o  the  filter  was  simulated  on  input  data  sets  already  contained  in 
.he  computer's  memory.  The  filter  used  two  inputs:  a  primary  input 
.hat  contained  noise  and  a  desired  signal  and  a  reference  input  that 
:ontained  noise  that  was  somehow  correlated  with  the  noise  in  the 
irimarv  input.  The  reference  input  was  filtere-"1  and  subtracted  from 
.he  primary  input  s-uppressing  the  noise.  The  filter  produced  signal  to 
toise  improvements  of  around  7*  dFJ  in  most  cases  and  worked  quite  well 
rith  and  adaptive  gain  of  5.^  •IT'’. 
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affect  of  background  noise  on  the  parameters  internal  to  the 
algorithm  and  this  is  considered  essential  to  basic  understanding 
of  the  problem. 
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cancellation  compared  to  noise  canceling  microphones  but  adaptive 
noise  cancellation  where  the  primary  input  is  a  noise  canceling 
microphone  v_s  a  single  noise  canceling  microphone.  In  related 
work  for  the  Air  Force  researchers  at  Bolt,  Beranek,  and  Newman, 
Inc.  [22]  have  experimented  with  various  multisensor 
conf igurations  and  have  found  a  combination  of  a  noise  canceling 
microphone  at  the  lips  and  an  accelerometer  attached  to  the 
throat  to  be  best.  The  accelerometer  is  insensitive  to 
environmental  noise  but  is  bandlimited  to  2  kHz.  The 
accelerometer  output  was  added  to  the  microphone  output  high  pass 
filtered  at  1200  Hz.  Also,  in  work  for  the  Navy,  Ketron,  Inc. 
[231  has  developed  a  second  order  gradient  noise  canceling 
microphone  (NC-104LF).  Their  tests  have  shown  substantially 
improved  performance  over  first  order  gradient  microphones. 

As  noted  in  the  introduction,  background  noise  is  a  more 
serious  problem  in  digital  voice  privacy  systems  employing 
the  LPC-10  speech  compression  algorithm  [12].  Therefore,  future 
efforts  should  incorporate  the  LPC-10  algorithm  and  measure 
intelligibility  after' synthesis .  Cadet  1/c  P.  A.  MENCEL  is 
presently  attempting  to  implement  LPC-10  in  FORTRAN  IV  on  an  LSI- 
11/2.  We  have  also  requested  the  PDP-11  FORTRAN  77  program 
referred  to  in  [12]  from  the  National  Security  Agency  and  we  hope 
to  modify  it  for  running  on  an  LSI-11.  These  programs  will  not 
run  in  real  time  as  required  for  DRT’s.  The  best  method  of  real 
time  implementation  would  be  a  TMS  320  based  coprocessor  board 
for  an  LSI-11.  An  alternate  method  would  be  to  acquire  two 
digital  voice  privacy  units  but  this  would  not  allow  access  the 
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weight  FIR  and  HR  filters,  plot  their  impulse  and  frequency 
responses  and  automatically  generate  TMS  320  assembly  language 
code  for  downloading  and  real  time  implementation.  The  program 
for  FIR  design  is  that  described  in  [21]  and  installed  on  DTSS 
when  PHOF  WOLCIN  joined  the  Academy  faculty  in  1984.  Cadet  1/c 
B.  M.  LAM  is  presently  modifying  the  program  to  produce  TMS  320 
code  for  real  time  implementation. 

V.  Conclusions  and  Recommendations  for  Future  Research 

We  have  demonstrated  that  under  some  conditions  digital 
filtering  can  improve  the  signal  to  noise  ratio  of  speech 
corrupted  by  background  engine  noise.  There  are  several 
questions  suggested  by  our  research  that  remain  unanswered. 

Does  improved  intelligibility  result  from  this  improved  SNR? 
Diagnostic  rhyme  tests  Cll]  are  necessary  to  resolve  this.  DRT’s 
are  quite  simple  and  straightforward  but  will  require  many  man¬ 
hours  of  effort.  Cadet  projects  may  be  the  best  method  of 
performing  DRT's.  How  intelligibility  is  measured  is  worth 
learning  and  the  data  analysis  is  a  good  exercise  in 
undergraduate  statistics. 

The  type  of  analysis  done  in  the  context  of  this  project  has 
not  lent  itself  to  a  good  comparison  of  adaptive  noise 
cancellation  v^  noise  canceling  microphones.  Again,  this  sort  of 
question  can  only  be  answered  using  DRT's.  The  general  consensus 
among  experts  in  the  field  is  that  noise  canceling  microphones 
are  essential  in  high  ambient  noise  environments  independent  what 
else  is  done.  Therefore  the  basic  question  is  not  adaptive  noise 


Conversations  with  other  researchers  in  the  field  indicate  that 
improved  signal  to  noise  ratio  may  not  necessarily  imply  improved 

intelligibility  and  that  DRT's  are  necessary  to  evaluate  a 

proposed  system.  The  human  auditory  system  is  very  sophisticated 
and  can  extract  information  from  highly  corrupted  signals  but 
non-linear  processing  may  destroy  cues  the  human  system  is  using 
not  improve  intelligibility. 

IT.  Publications  and  Presentations 

Presentations  of  papers  based  the  research  have  been  made  at 
two  IEEE  sponsored  conferences  and  one  local  IEEE  meeting  to  date 

and  a  fourth  abstract  has  been  submitted  for  presentation  at  a 

future  conference.  As  noted  in  [16]  M.  D.  SAKAHARA  presented  his 
paper  "Adaptive  Tioise  Cancellation"  at  Electro  '84'  at  Boston  in 
May  1984.  He  also  presented  "Real  Time  Adaptive  Noise 
Cancellation"  at  the  April  1984  IEEE  New  London  Subsection 
meeting.  This  paper  was  included  in  [16]. 

B.  B.  PETERSON  presented  "Audio  Frequency  Adaptive  Filtering 
Using  the  TMS  320  Microprocessor"  at  the  IEEE  Digital  Signal 
Processing  Workshop  at  Chatham,  MA  in  October  1984.  In  addition 
K.  U.  DYKSTRA  and  M.  D.  SAKAHARA  were  coauthors  of  the  published 
summary  and  LCDR  DYKSTRA  attended. 

An  abstract  titled  "An  Integrated  Approach  to  the  Design  and 
Implementation  of  Digital  Filters"  has  been  submitted  to  the 
IEEE/ASEE  Frontiers  in  Education  Conference  at  Golden,  CO  in 
October  1985.  In  addition  to  the  BASIC  program  in  Appendix  C, 
the  paper  will  describe  interactive  programs  to  design  fixed 
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filtering  the  reference  input  with  an  80  weight  filter  and 
subtracting.  The  sampling  frequency  was  8  kHz. 

The  program  was  then  modified  to  sample  a  third  channel  which 
was  the  audio  output  of  an  FM  tuner.  In  software  this  signal  was 
then  amplified  and  added  to  the  primary  input  and  resultant 
signals  run  through  the  adaptive  filter  algorithm.  Appendix  F  is 
this  FORTRAN  program.  A  listing  of  the  graphics  package  written 
to  drive  the  DEC  Gigi  Graphics  Terminal  and  the  HP  7470A  Plotter 
is  not  included  because  of  it's  length  and  because  the  subroutine 

I 

calls  are  obvious.  This  package  was  the  most  extensive 
programming  effort  undertaken  in  the  context  of  this  research  and 
the  excellent  work  of  LTJG  W.  M.  DUPRIEST-is  acknowledged. 

Unlike  previous  efforts,  it  is  now  possible  to  know  exactly 
what  is  signal  and  what  is  noise.  Figure  2  is  a  sample  plot  of 
the  results.  The  signal  to  noise  ratio  has  been  improved 
from  -3-3  db  to  +6.0  db.  Improvements  in  signal  to  noise  ratio 
were  consistently  near  10  db.  Reductions  in  noise  power  of  20  db 
as  noted  when  no  signal  was  added  (Figure  1)  are  not  expected 
because  what  is  signal  to  the  overall  filter  algorithm  is 
actually  noise  to  the  weight  update  algorithm.  The  fluctuation 
in  weights  caused  by  the  larger  output  degrades  the  filter 
performance.  This  effect  is  minimized  by  reducing  the  adaptive 
gain  at  the  expense  slower  convergence  and  poorer  tracking  in 
time  varying  conditions. 

What  these  recent  results  do  indicate  is  that  adaptive 
filtering  can  work  very  well  when  the  noise  is  highly  correlated. 


e.  Frequency  Domain  Adaptive  Filtering 

Just  as  their  fixed  weight  counterparts ,  frequency  domain 
adaptive  filters  require  only  a  fraction  of  the  arithmetic 
operations  required  in  a  time  domain  filter  of  equivalent 
complexity.  The  frequency  domain  adaptive  filter  algorithm 
described  in  [5]  was  implemented  in  FORTRAN  for  filter  lengths  of 
up  to  512  samples.  Due  to  filter  length  and  good  convergence 
properties  the  filter  performance  was  better  than  that  of  the 
adaptive  filters  in  a.  and  c.  above  but  not  as  good  as  the  notch 
filter  in  d.  Detailed  analysis  of  the  filter  is  contained  in 
Appendix  E. 

DRT's  using  a  real  time  version  of  the  algorithm  would  be 
necessary  to  make  a  more  complete  evaluation.  The  TMS  320  is 
fast  enough  but  due  to  limited  memory,  implementation  wifch  our 
existing  hardware  is  not  possible. 

f.  Adaptive  Filtering  of  Automobile  Engine  Noise 

Our  most  recent  effort  was  to  make  a  two  microphone  recording 
of  engine  noise  from  a  Datsun  710  with  no  muffler  running  at  idle 
(approximately  1000  rpra)  in  the  Power  Engineering  Laboratory. 

Due  to  the  faulty  exhaust  system,  it  was  necessary  to  operate  the 
ventilation  system  at  full  power  adding  another  significant 
component  of  noise. 

When  the  data  was  plotted  in  the  laboratory,  the  two  signals 
were  seen  to  be  much  more  correlated  than  had  been  observed  in 
previous  recordings.  Figure  1  shows  two  signals  and  the  results 
of  adaptive  filtering  with  no  voice  signal  added.  98.9%  of  the 
power  in  the  primary  input  has  been  canceled  by  adaptively 


Results  and  Conclusions 


?ura  signals  (  sins,  square  and  triangle  )  wer»  used  to 
simulate  noisa  and  desired  signal.  A.  square  wave  was  used  to 
simulate  noise  sines  it  contained  the  most  harmonics  and  could  be 
shaped  to  any  of  the  noise  signals  placed  in  the  primary  input. 

Three  ',/avetek  signal  generators  were  used  to  supply  the 
noisa,  signal,  and  clock  pulses  for  the  sampler.  Hewlett-Packard 
frequency  counters  and  oscilloscopes  were  used  to  monitor  the 
inputs  and  outputs  of  the  computer  simulation.  In  the 
simulations,  the  noise  in  the  reference  and  primary  inouts  './ere 
of  the  same  frequency  and  phase.  This  simplifies  predicting  the 
impulse  response  the  filter  should  converge  to  allowing  for  an 
evaluation  of  the  filter's  performance. 

The  primary  input  contained  a  square  wave  as  the  desired 
signal.  This  enables  the  calculation  of  signal  to  noise  ratios  ( 
SUP.  )  based  on  the  amplitudes  of  the  noise  and  the  desired 
signal.  The  SJJR  improvement  of  the  filter  could  be  estimated  and 
the  performance  of  the  filter  evaluated  based  on  the  31JR 
improvement . 


A  triangle  wave  represented  the  noise  in  the  primary  input. 
This  simplified  predicting  the  imp'ulse  response  the  filter  should 
converge  to.  The  Fourier  series  of  -a  square  of  unity  amplitude 
is  .  ..  .  * 

4/V  (  sin  wt  +  1/3  sin  3wt  +  1/5  sin  5wt  +  .  .  .  ) 

and  a  triangle  wave  of  same  amplitude 

3/r  (  sin  wt  -  1/9  sin  3wt  +•  1/25  sin  5wt  -...). 

In  order  to  change  the  square  wave  in  the  reference  input  into 
the  triangle  wave  in  the  primary  input,  the  square  wave  must  be 
multiplied  by  another  square  wave  of  the  same  frequency. 

‘/hen  the  weights  were  plotted,  this  indeed  seemed  to  be  the 
result.  There  were  small  deviations  in  the  imoulse  response  of 
the  filter  and  was  evident  in  the  incomplete  cancellation  of  the 
noi3e  in  the  primary  input.  Although  noise  cancellation  was.  not 
total,  5HR  improvements  of  up  to  32  dB  were  observed. 

The  adaptive  coefficient  was  varied  and  had  a  proportional 
affect  on  the  convergence  time  of  the  filter,  '/hen  u  was  1.25 -IT" 
,  the  filter  converged  in  12  msec,  '/hen  u  was  4.d  •  l^"1  ,  the 
convergence  time  was  45  msec  or  it  took  about  four  times  as  long. 
In  both  cases,  SUP.  improvements  of  2n  dB  was  observed. 


Then  u  was  increased  beyond  2 .  h  •  lh'*,  the  filter  seemed  to 
have  no  effect  on  the  orinary  input.  Since  the  adaptive 
coefficient  was  so  large,  the  -./eights  constantly  over  shot  the 
desired  solution  and  did  not  reduce  the  noise  levels  in  the 
primary  input.  The  adaptive  gain  was  not  sufficiently  high 
enough  to  cause  instability,  ju3t  large  enough  to  cause  excessive 
oscillations.  ’Then  the  adaptive  gain  was  increased  beyond  2.’, 
the  filter  quickly  became  unstable  and  blew-up.  This  was  the 
threshold  where  the  weights  were  allowed  swing  so  far  that  they 
overcame  the  reducing  tendencies  of  u. 

The  advantages  of  an  adaptive  filter  were  especially  evident 
in  the  case  whan  the  noise  and  signal  were  similar  in  frequency. 

A  run  -with  the  noise  at  5h'1'  Hz  and  the  desired  signal  at  7  5h  Hz 
produced  a  SNR  improvement  of  32  13  i 

The  LMS  filter  seemed  to  suffer  when  the  noise  was  at  a 
significantly  lower  frequency  than  the  signal.  The  filter  would 
glitch  "  every  time  the  triangle  wave  would  change  slope.  The 
filter  would  quickly  converge  to  the  desired  signal  while  the 
slope  was  constant,  however  would  loose  track  when  the  slope 
changed.  this  is  due  ta  the  fact  that  the  LMS  filter  follows  the 
gradient  of  the  noise.  A  low  frequency  triangle  wave  has  a  low 
gradient  while  scribing  constant  slope,  however  the  gradient 
greatly  increases  when  the  triangle  wave  changes  slope. 

Over  all,  the  LMS  algorithm  performed  quite  wall  in  this 
application.  It  produced  SNR  improvements  of  around  2hd3  and 
exceeded  30  dS  in  some  cases.  The  adaptive,  coefficient’  was 
critical  in. the  performance  of  the  system.  The  lower  frequencies 
tended  to  require  higher  adaptive  gains  to  optimize  the  filter, 
while  the  higher  frequencies  required  lower  adaptive  gains.  The 
adaptive  gains  depended  mostly  on  the  contents  of  the  reference 
input  and  the  noise  frequencies.  The  adaptive  gains  seemed  to  be 
independent  of  the  desired  signal.  An  adaptive  gain  of  S.'VlT4 
•worked  the  best  for  all  cases  of  noise  and  signal  combinations. 
The  LMS  filter  worked  well  for  most  combinations  of  noise  and 
desired  signals. 


Design  of  interface  between  L3I-11/2  and  TRW  s  TDC-ldldj 


A  cookbook  approach  was  used  in  designing  the  interface 
between  DR’?*  s  TDC-ihinj  Multiplier  Accumulator  chin  and  the  LSI- 
11/2  computer.  MDS  Systems  general  purpose  interface  module  ( 
MS:-17n  )  was  used  to  mount  the  TDC-ldlGJ  and  supporting 
hardware . 

The  interface  module  provided  input, output  and  address  ports 
that  were  common  on  the  Q— bus  of  the  LSI— 11.  This  simplified  the 
decoding  and  interface  logic  greatly  since  bus  protocol  did  not 
have  to  be  considered.  It  was  determined  that  the  TDC-l^l^J 
could  multiply  and  accumulate  with  in  the  tine  of  a  DATO  cycle  of 
the  LSI-11.  Thus  the  only  timing  consdierations  made  were  in 
delaying  signals  going  into  the  TDC-171'’. 

The  Y  multiplicand  and  the  least  significant  IS  bits  of  the 
output  share  pins  on  the  TDC-l^lD.  Logic  had  to  be  provided  to 
isolate  the  input  and  output  ports  on  the  interface  module. 
‘.Ion inverting  tri-state  bus  transceivers  (  Ti  ‘  s  74LS243  )  were 
used  to  isolate  the  ports.  The  3  side  of  the  transceivers  were 
connected  to  the  TDC-lDldJ  and  the  A  side  to  the  respective  I/O 
location.  The  control  inputs  for  the  transceivers  are  Gab  (L) 
and  Gba  (El).  EThen  Gab  is  asserted,  the  data  is  allowed  to  flow 
from  the  A  side  to  the  3  side.  When  Gba  is  asserted,  the  data  is 
alowed  to  flow  from  the  B  side,  to  the  A  side.  Asserting  both 
control  inputs  could  result  in  destructive  oscillations  because 
the  transceiver  is  trying  to  pass  data  in  both  directions  at  the 
same ‘time.  Mot  asserting  both  control  incuts  isolates  side  A 
from  side  3. 

In  the  case  of  a  Y  input,  the  Gba  input  is  held  low  not 
asserting  it  and  Gab  was  asserted  when  the  input  data  was  valid. 
The  reason  Gba  was  held  low  was  to  prevent  the  possibility  of 
sanding  the  transceiver  into  destructive  oscillations.  By 
keeping  Gba  low,  the  transceiver  can  be  alternated  from  isolation 
to  passing  data  from  A  to  3  (  the  desired  direction  ) .  The  logic 
for  the  Tab  input  was  (  ODOUT  •  ADD3  )  (L).  ADD3  is  the  address 
location  indicating  a  Y-input  on  the  bus  and  BDOUT  is  asserted 
when  the  data  in  the  input  port  is  valid.  The  Tab  signal  was  used 
to  supply  the  CLX  Y  signal.  The  three  inverters  delay  the  CLX  Y 
signal  to  account  for  propagation  delay  of  the  transceivers. 

Since  the  X-input  shares  the  input  port  with  the  Y-input  and 
control  byte,  it  was  isolated  from  the  input  port  when  not  in 
use.  This  was  accomplished  by  using  tri-state  noninverting  bus 
drivers  (  Ti’s  74367  ).  This  circuitry  i3  net  absolutely 
necessary  because  the  X-input  is  Loaded  on  to  the  TDC-L'HGJ  only 
when  CLK  X  is  asserted.  But  was  added  as  a  precautionary  measure 


only.  The  control  logic  for  the  bus  drivers  is  (  ADD'’  *  3D0UT  ) 
(L).  ADD  2  is  the  address  location  indicating  an  X-input  will 
appear  on  the  input  port  and  3DCUT  is  asserted  when  the  data  on 
the  input  port  is  valid.  This  control  signal  was  U3ad  to  supoly 
the  CLK  X  signal  with  inverters  to  delay  to  3ignal  to  account  for 
propagation  delays  in  the  bus  iriver3. 

The  7DC-171DJ  has  five  data  control  bits: 

Accumulate  (  ACC  )  bit  00 
Subtract  {  SUB  )  bit  "*1 
Preload  (  PRSL  )  bit  02 
Two's  Complement  (  TC  )  bit  03 
Bound  off  (  BUD  )  bit  ^4 

all  asserted  high.  It  was  decided  for  versatility  that  these 
control  bits  be  controlled  by  the  LSI-11.  The  control  bits  are 
loaded  onto  the  TDC-1010J  at  different  times  makinq  it  difficult 
or  impossible  to  send  them  with  the  data  going  into  the  board. 
Type  0  flip  flops  (  Ti's  7474  )  are  used  to  store  the  control 
byte.  The  clock  signal  for  the  flip  flops  is  (  ADD1  •  9 DOLT  )  (ft). 
ADD1  is  the  address  location  indicating  that  the  control  byte  is 
to  appear  on  the  input  port  and  3DOrjT  is  asserted  when  the  data 
on  the  input  port  is  valid. 

The  output  of  the  TDC-lHldj  is  divided  into  three  16  bit 
words:  the  least  significant  product  (  L3P  ),  most  significant 
product  (  MSP  )  and  extended  control  product  (  XTP  ) .  The  LSP  is 
time  shared  wit  the  Y-input  and  similar  circuitry  is  used  to 
isolate  the  input  and  output  ports  as  in  the  Y-input.  The  tri¬ 
state  output  of  the  TDC-l'31'’iJ  is  controlled  by  throe  inputs  TSX, 
TSM  and  TSL.  PP.EL  must  be  held  low  if  the  TDC-1710J  is  to 
canerate  output.  The  control  inputs  are  asserted  low  and  only  one 
of  them  may  be  asserted  at  a  tine  for  this  design.  A  separate 
address  has  been  allocated  for  each  output  word  with  TSM  clocking 
the  output  to  the  output  registers.  The  output  of  the  address 
decoder  on  the  interface  module  is  asserted  low  and  is  connected 
directly  to  the  appropriate  output  control  pins.  The  interface 
module  requires  an  input  enable  signal  to  indicate  that  the  data 
on  the  output  port  is  valid.  This  is  accomplished  by  ORing  the 
output  control  signals  and  using  an  inverter  to  account  for 
propagation  delays. 

Six  address  locations  were  used  to  implement  the  TDC-l^l^J  : 


Address 

"unction 

number 

Location 

7 

CTRL  2yta 

17  'I'jnn 

X-  Input 

2 

Y-  Input 

177^72 

3 

LSP  P rod. 

17^712 

4 

MSP  Prod. 

17',774 

5 

XTP  Prod. 

177714 

The  unusual  order  of  address  locations  is  due  to  the  address 
decoder  on  the  interface  nodule.  The  decoder  an  A17D  gate  for  3 
significant  bits  and  a  BCD  to  Decimal  (  Ti ’ s  7442  )  for  the  lower 
three  bits.  3its  '"l  and  72  are  hard  wired  to  the  3  and  C  inputs 
and  the  output  of  the  AiTD  gate  is  wired  to  the  D  input  of  the 
decoder.  3it  73  was  wired  to  the  A  input  and  resulted  in  the 
unusual  address  order.  3it  73  was  brought  to  the  decoder  to 
enable  it  to  decode  eight  address  instead  of  four. 

To  use  the  TDC-1717J  ,  the  output  registers  must  be  zeroed. 
This  can  be  accomplished  by  sending  a  7  to  the  X- input  and 
multiplying  without  accumulation.  The  MS?  must  be  outputted  to  a 
dummy  location  to  clock  the  output  registers.  The  other 
functions  of  the  TDC-1717J  can  then  be  executed  by  altering  the 
control  byte  as  necessary.  'Then  ever  output,  accumulation  or 
subtraction  is  desired,  the  MSP  must  be  moved  to  a  dummy  location 
since  .asserting  the  MSP  ‘clocks  the  output  clock  and  activates  the, 
accumulator  and  clocks  the  output  registers. 
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Appendix  I.  Flow  Chart  and  Coded  Program 
for  LMS  Adaptive  Filter. 


FORTRAN  IV 


'.'02 . 5 


rut  os-Dec-B3  oo:oo:oo 


°AGE  001 


0001 


0002 

0003 


000-4 

0005 

0006 

0007 

0008 

0003 


0010 


0011 

0012 

0013 

0014 


001? 

0016 

0017 

0013 

0019 

0020 

0021 


0022 


C  PROGRAMMER:  Michael  D.  SaKahara 
C 

C  DATE!  November  13.  1963 
C 

PROGRAM  ADAPT 

C  THIS  PROGRAM  IMPLEMENTS  AN  ADAPTIVE  ,rILTER  ALGGRYTHM.  THE  FILLER 
C  IS  IMPLEMENTED  IN  THE  BEGIN ING  OF  THE  PROGRAM  AND  THE  USER  INPU” 

C  PORTION  OF  THE  PROGRAM  APPEARS  AT  THE  END  GF  THIS  PROGRAM. 

C 

C 

C  INSTRUCTIONS:  WHEN  ANSWERING  YES  TO  ANY  OF  t-jc  QUESTIONS.  ENTER 
C  ANY  VALID  INTEGER .  HOWEVER  WHEN  ENTERING  NO  ENTER  A  '1'  ONLY. 

C 

c 

c 

DIMENSION  WEIGHT ( 100 ) »  INI (3000).  IN2(?000>.  REAL (201 ) ,  °“A(201) 
GOTO  50 
C 

C  ADD  DELAY  TO  AVOID  INTERRUPT  FROM  KEYBOARD 
C 

10  DO  100  I » 1 r  10000 
100  CONTINUE 

DO  110  1*1.5000 
INK  I  )*0 
IN2  < I ) *0 
110  CONTINUE 
C 

C  GET  NUMS  SAMPLES  FROM  A/D  CONVERTER 
C 

CALL  XTSMP2  ( NUMS. NC 1 .NC2. INI , IN2) 

C 

C  ADD  DELAY  TO  AVOID  INTERRUPT  FROM  KEYBOARD 
C 

20  DO  120  1*1,  10000 
120  CONTINUE 

DO  130  I*  NUMW+1 .NUMS-NUMD 
SUM*0 . 0 

c 

C  IMPLEMENT  FIXED  WEIGHT  FILTER 
C 

DQ  140  J-l.NUMW 

sum*sum+we:ght( j  >  *  in 1 ( :-j  > 

140  CONTINUE 

K*IN2< I+NUMD) 

OUT  *K-SUM 

WRITE  (7,800)  K , CUT 
300  FQRMAT  ( IX , I8.G12 . 5 ) 

C 

C  OUTPUT  OF  FIXED  WEIGHT  FILTER  TO  D/A  (INPUT,  CUTOUT) 

C 

CALL  DAOUT (K,INT(QU7) ) 

C 

C  IMPLEMENT  ADAPTIVE  FILTER 
C 


A  15 


®OR"RAN  IV 

V<'»2.5  Tue  OS-Dec-83  00 1 00 100  aAGE  002 

0023 

<-:sADAPT*CUT 

002* 

DC  1 50  J  =  i  ,  NUMU 

0025 

WEIGHT  (  J )  *  WEIGH"  (J)  ■*■(.-••*  I N  i  ■'  I  — -  .■  > 

\V, 

0025 

1 50 

CONTINUE 

<•'027 

120 

CONTINUE 

0029 

CAi_U  DAOUT  ( 0  ,  0  > 

~  '' 

002S 

READ  ( 5 ,  905 )  IANS 

0030 

r> 

IF  (  IANS  .EG.  0  )  GO"Q  25 

L 

C  aL 

r» 

.OT  WEIGHTS  ON  STRIP  CHAR" 

0032 

U 

DO  ISO  I*1,NUMW 

0033 

IF  <A8S  (WEIGHT  (I))  . uT .  WHAM  )  WMAX*  ASS  (WEIGH"  (I)) 

• 

0035 

ISO 

CONTINUE 

0035 

WRITE  (7,810)  WMAX 

0037 

WSCL*500 .0/WP1AX 

*  ' 

0038 

DO  155  1*1,  NUMW 

0039 

WRITE  (7,800)  I,  WEIGHT  (I) 

0040 

TEMP*WEIGHT  (I)  *  WSCL  +  0.5 

00*1 

CALL  DAGUT  ( I NT  (TEMP),  0) 

00-2 

1S5 

CONTINUE 

N* 

00*13 

CAUL  DAOUT  (0,0) 

v-l’ 

L 

-  C  PL 

,QT  FREQUENCY  RESPONSE  Or  FILTER  (MAGNITUDE  AND  °HASE> 

• 

0644 

c 

NF*2Q0  *  • 

0045 

CAUL  FREGMP  (WEIGHT,  NUHW,  NF,  REAL,  °HA ) 

0045 

RMAX*0 . 0 

■  \ 

0047 

®MAX*0 . 0 

0  0*8 

RMIN* 10000 . 0 

0043 

PM  IN* 10000.0 

**4 

0050 

DO  170  I*  1,  NF 

* 

0051 

IF  (  REAL ( I )  .GT.  RMAX  )  RMAX*REAL ( I ) 

'  1 

0053 

IF  (  REAL ( I )  . LT .  RMIN  )  RMIN*REAL ( I > 

0055 

IF  <  P«A(I)  .GT.  PMAX  )  °MAX*°HA ( I ) 

0057 

IF  (  PHA(I)  . L" .  PM IN  )  PMIN*  °HA ( I > 

0059 

170 

CONTINUE 

0050 

RSCL=RMIN*(-i .0) 

0051 

I F  (  RMAX  .GT.  C-:.0>*RMIN  )  RSCL*RMAX 

<5053 

°SCL*PMIN* ( -1 .0) 

005* 

IF  (  PMAX  . GT .  (-1.0)*®“ IN  )  ®SCL  =  0MAX 

0055 

®SCL*500 . 0/°SCL 

0057 

RSCL*500 . O/RSCL 

0  053 

DO  180  1*1,  NF 

0059 

REAL ( I ) * ( REAL ( I ) ♦RSCL ) *0 . 5 

i 

0070 

N1*INT(REAL( I ) ) 

*„  • ■„ 

0071 

PHA ( I ) * ( ®UA( I >*?SCL)*0.5 

0072 

N2* INT ( PHA ( I ) ) 

.  -  A 
■  1 

0073 

WRITS  (7,810)  REAL ( I ) , PHA ( I ) 

007* 

810 

FORMAT  ( 1X,G12.5, 1X/G12.5) 

0075 

CALw  DAOUT  (N1,N2) 

_  ' 

0075 

180 

continue 

0077 

CALL  DAOUT  (0,  0) 

*.*/ 

c 


GR7RAN  IV 


VO  2 


<ie  QS-Dec-83  OOIOOIOO 


C  SET  U®  -QP  RERUN  I~  DESIRED 
C 

079  WRITE  (7,900) 

073  SCO  FORMAT  (IX,  'ANOTHER  GO.*  1»NC  ?  ‘  ,  S  > 

050  READ  (5,205)  IANS 

OSl  305  -ORmaT  Cl) 

032  Ic  (IANS  .EG.  1)  GGTQ  1000 

OS*1  25  WRITE  (7,S10) 

035  9  ■  0  FORMAT  (IX,  'GET  NEW  INPUTS  ?',$) 

03G  READ  (5,305)  IANS 

087  WRIT?  <7,  SJ5) 

088  915  FORMAT  (IX,  'RESET  WEIGHT  VECTOR  ?',?) 

083  READ  (5,305)  IANS1 

090  IF  (IANS1  .EG.  1)  GOTO  GO 

092  DC  ISO  I  *  1,  NUMW 

092  WEIGHT ( I ) *0 

094  130  CONTINUE 

095  30  WRITE  (7,320) 

0S3  920  CGR!"*AT  (1>(,  'NEW  ADAPTIVE  CO-EF- ICSIN'  ?',«) 

097  READ  (5,  905)  I ANSI 

098  IF  (IANS1  .EG.  1)  GOTO  40 

100  WRITE  (7,  925) 

101  925  FORMAT  (IX,  'INPUT  NEW  ADAPTIVE  WEIGHT  CGEF^ICIE! 

102  READ  (5,930)  ADAPT 

103  930  FORMAT  (G14.7) 

104  a. o  IF  (IANS  .EG.  1)  GOTO  20 

106  GOTO  10 

C 

C  ACTUAL  3EGINNING  OF  PROGRAM  (  FIRST  RUN  ) 

C 


39T9AN  IV  Storaae  Ma®  for  srosram  unit  ADA?" 


ocai 

Var i as  I 

«S,  .  PScC” 

5DAT4  r 

Size 

=  073174  ( i 

.052 .  K*t c 

r  d  5  ) 

awe 

T':*P  9 

Offset 

Name 

T  y  pe 

Offset 

'■ame 

T  y  p  e 

Offset 

DA?T 

9*4 

053052 

LI 

9*4 

073046 

I  *2 

073016 

ANS 

1*2 

073076 

I  AMS  1 

1*2 

073132 

J 

1*2 

073036 

1*2 

073040 

NCI 

I  *2 

073022 

NC2 

1*2 

073024 

r 

1*2 

073074 

NUMD 

—  "  tmm 

073030 

MUMS 

1*2 

053020 

UMM 

1*2 

07302S 

N 1 

1*2 

073126 

M2 

1*2 

073130 

LT 

<9*4 

073042 

°MAX 

9*4 

073102 

l3M  I  N 

9*4 

053 112 

SCL 

9*4 

073 1 22 

9MAX 

9*4 

073076 

9MIN 

9*4 

073106 

SCL 

9*4 

073116 

SUM 

9*4 

073032 

TgMp 

9*4 

073070 

MAX 

9*4 

073060 

WSCL 

9*4 

073064 

oca  I 

and  COMMON  Arrays 

• 

• 

ame 

T  ype 

1  Section 

Offset 

_ 1 _ 

~'S  i  z  e - 

Dimens ions 

N 1 

1*2 

*DATA 

000620 

0234 

20  <  7000 . ) 

( 5000 ) 

N2 

1*2 

SDATA 

024240 

0234 

20  (  7000 . ) 

( 5000  > 

HA 

9*4 

SDATA 

071324 

00144a  <  402.) 

(201  > 

SAL 

9*4 

SDATA 

047660 

001444  (  402.) 

(201 

SIGHT 

R*4 

SDATA 

000000 

0006 

20  (  200 . ) 

(  100) 

uSro'JtineSf 

fund  ions  » 

Statement  and  Processor 

-Defined 

fund 

ions! 

awe 

Type 

Name  •  Ty 

pe  Name  T 

ype  .Name  . 

Ty*« 

Name 

T  y»j 

SS 

9*4 

DAOUT  R*4  r 9EQMP 

9*4  INT 

1*2 

XTSMP 

2  9*4 

A18 


VARIABLE  LIST 


ADAPT-  Adaptive  coefficient 

H  -  Product  of  the  filter  output  and  adaptive  coefficient 

I  -  Pointer  for  Do  loops 

IANS  -  Variable  for  responses  from  the  keyboard 
IANS1-  Variable  for  responses  from  the  keyboard 
J  -  Pointer  for  Do  loops 

K  -  The  present  input  value  with  delay  accounted  for 

NCI  -  Address  pointer  for  the  channel  of  the  noise  source 
NC2  -  Address  pointer  for  the  channel  of  the  signal  +  noise  source 
NF  -  Number  of  frequencies  between  0  and  fs/2  the  frequency  spectrum 
of  the  weight  vector  is  calculated 
NUMD  -  Number  of  delays  in  the  filter 
NUMS  -  Number  of  samples  taken 
NUMW  -  Number  of  weights  in  the  filter 

N1  -  Integer  value  of  array  REAL  scaled  between  0  and  500 
N2  -  Integer  value  of  array  PHA  scaled  between  0  and  500 
OUT  -  Filter  output 

PMAX  -  Maximum  value  of  the  array  PHA 
PMIN  -  Minimum  value  of  the  array  PHA 

PSCL  -  Factor  required  to  scale  array  PHA  to  range  from  0  to  500 
RMAX  -  Maximum  value  of  the  array  REAL 
RMIN  -  Minimum  value  of  the  array  REAL 

RSCL  -  Factor  required  to  scale  the  array  REAL  to  range  from  0  to  500 
SUM  -  Summation  variable  for  the  fixed  filter 


TEMP  -  Scaled  value  of  the  array  WEIGHT 
WMAX  -  Maximum  value  of  the  array  WEIGHT 

WSCL  -  Factor  required  to  scale  the  array  WEIGHT  to  range  from  0  to  500 


ARRAY  LIST 

INI  -  Input  array  containing  noise  samples  (5000) 

IN2  -  Input  array  containing  signal  +  noise  samples  (5000) 

PHA  -  Phase  values  of  the  frequency  spectrum  of  the  weight  vector  (201) 
REAL  -  Magnitude  values  of  the  frequency  spectrum  of  the  weight  array (201) 
WEIGHT-  Weight  vector  (  100  ) 
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W(k+1 )  =  W.(k)  +  uYS(k-i) 

L  L 


where  W.  are  the  N  weights,  u  is  the  adaptive  gain  coefficient,  Y 
is  the  last  output  of  the  system,  and  S  are  the  reference  inputs 
(Na).  The  sequence  of  four  instructions  that  performs  this 
operation  are: 


ZALH  * 

MPY  63 
APAC 
SACH  »- 
ZALH  * 

MPY  67 
APAC 
SACH  *- 
ZALH  * 

MPY  66 
APAC 
SACH  »- 
etc . 

The  ZALH  *  instruction  zeroes  the  accumulator  and  loads  the 
contents  of  the  data  memory  location  pointed  to  by  the' current 
auxiliary  register  (set  to  point  to  the  weights)  into  the  upper 
16  bits  of  the  accumulator.  The  MPY  multiplies  the  contents  of 
the  data  memory  location  (location  63,  67  ...  contains  the  past 
reference  inputs)  by  the  T  register  (which  contains  the  adaptive 
gain  coefficient  and  output  product  uxY).  The  APAC  instruction 
adds  the  contents  of  the  P  register  (result  of  last  multiply)  to 
the  accumulator  and  stores  result  in  the  accumulator.  The  SACH  *- 
instruction  stores  the  upper  16  bits  of  the  accumulator  into  the 
data  memory  location  pointed  to  by  the  current  auxiliary  register 
(still  points  to  the  location  of  the  weight)  and  the  auxiliary 
register  is  decremented. 


C3 

The  tapped  delay  line  portion  of  the  program  calculates: 

M 

F(k)=  W.(  k ) S ( k- i ) 

i  =  1  *• 

where  F  is  the  output  of  the  adaptive  filter,  N  is  the  number  of 
weights,  S  are  the  past  N  reference  inputs,  and  W  are  the  M 
weights.  This  can  be  implemented  with  a  sequence  of  two 
instructions . 

LTD  67 
MPY  »- 

LTD  66  1 

MPY  *- 

LTD  65 

MPY  »- 

etc . 

The  LTD  instruction  places  the  contents  of  the  data  memory 
location  (locations  67,  66,  ...  contain  the  past  reference' 
inputs)  in  the  T  register  to  set  up  for  the  next  multiply,  add 
the  result  of  the  last  multiply  which  is  in  the  product  register 
to  the  accumulator,  and  shifts  the  data  memory  to  the  next  data 
memory  locationC  contents  of  location  67  would  be  placed  in  63). 
The  MPY  *-  will  multiply  the  contents  of  the  T  register  by  the 
contents  of  the  data  memory  that  the  current  auxiliary  register 
points  to  (it  is  set  up  to  point  to  the  data  memory  containing 
the  weights)  and  then  decrement  the  auxiliary  register  to  point 
to  the  next  weight. 

The  adaptive  filter  portion  of  the  program  uses  a  series  of 
four  instructions.  The  new  weights  are  calculated  by: 


C2 


Figure  2  shows  a  flow  diagram  of  the  a  TMS-320  program  that 
performs  the  adaptive  noise  cancellation.  The  program  begins  with 
an  initialization  portion  to  set  up  certain  constants.  In  the 
main  loop,  there  are  several  steps  performed.  First  a  reference 
input  (N0)  is  read.  Next  the  tapped  delay  line  filter  is 
performed.  The  output  of  the  system  is  next  computed  and 
transferred  to  the  output  port.  Then  the  signal  plus  noise  input 
(S  +N})is  read.  Finally,  the  filter  weights  are  adaptively 
adjusted.  This  loop  is  performed  once  every  sampling  period. 

Delay  loops  are  also  part  of  the  loop  to  control  the  sampling 
rate.  To  show  the  power  of  the  TMS-320  instruction  set,  the 
adaptive  filter  portions  of  the  program  are  explained. 


Figure  2.  Flow  diagram  of  adaptive  noise  cancellation 
program. 


Appendix  C. 


Implementation  of  an  Adaptive  Noise 
Cancellation  Algorithm  on  the  TMS-320 


An  adaptive  noise  cancellation  system  has  been  developed  on 

the  Texas  Instruments  TMS-320.  A  block  diagram  of  the  system  is 

shown  in  Figure  Cl.  The  TMS-320  can  implement  a  68  weight 

adaptive  filter  operating  at  a  sampling  frequency  of  10.7  kHz. 

The  TMS-320  has  an  instruction  set  tailored  to  do  digital 

signal  processing.  With  its  200  nsec  instruction  cycle,  it 

becomes  a  very  powerful  processor  for  applications  such  as  this. 

The  program  to  do  the  adaptive  noise  cancellation  only  requires 

six  instruction  per  weight  or  1.2  usee  plus  the  overhead  needed 

to  do  input,  output,  and  minor  data  manipulations.  The  maximum 

% 

number  of  filter  weights  that  th'e  TMS-320  can  support  without 
major  hardware  modification  is  68. 


•  '.'I 

■*  ‘v*1 

Figure  Cl.  Adaptive  Noise  Cancellation  System.  .  ' 
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Enclosure  3  to  Appendix  D  Assembly  language  sampling  subroutine  called  fran 
rain  program.  (Also  used  with  programs  in  Appendixes  E  and  ?.) 

.TITLE  XTSMP2 
.  GL08L  XTSMP2 

f 

’  8.  B.  PETERSON  17  NOV  83 

J  SUBROUTINE  CALLED  BY  "CALL  XTSMP2 ( N , N 1 , N2 , IX 1 , I X2 ) " .  N=  #  OF 

J  SAMPLES  IN  EACH  OF  CHANNELS  NCI  AND  NC2  ,  N 1  *25S*nC  1  l  S  ETC.) 

;  AND  1X1  AND  1X2  ARE  THE  RETURNED  DATA  VECTORS. 

J  MAXIMUM  TRIGGER  RATE  IS  21,000  SAMPLES/SEC  (10,500  PER  CHANNEL 

• 

r 

IOC* 170400  ; ADDRESS  OF  A/D  CONTROL  REGISTER 

101*170402  ' ADDRESS  OF  A/D  OUTPUT  REGISTER 

XTSMP2 I  TST  i R5 ) 

MOV  0<R5)+,RO  JRO  *  #  OF  SAMPLES 
MOV  3<R5)«-,R1  !R1  *  N1 

MOV  0<R5)+,R2  JR2  *  N2 

MOV  ( R5 )  *• » R3  JR3  *  STARTING  ADDRESS  OF  DATA  VECTOR  1X1 

MOV  (R5)+,R4  ;R4  a  STARTING  ADDRESS  OF  1X2 

MOV  R1,@#I0C  *, ENABLE  EXTERNAL  TRIGGER,  SET  ADD  TO  CH  NCI 

LOOP:  BIT  #2OO,0#IOC  ?TEST  DONE  BIT 

SEQ  LOOP 

MOV  8#I0I,R5  " RESET  DONE  SIT,  THROW  AWAY  DATA 

TEST  1 :  BIT  #2OO,0#IOC  JTEST  DONE  BIT 

BEG  TEST 1  ;WAIT  AND  TEST  AGAIN  IF  NOT  SET 

MOV  R2,0#IOC  ;SET  ADDRESS  TO  CH  NCZ 
MOV  0#IOI , vR3 )+  READ  A/D  AND  PUT  IN  MEMORY 
LQ0P2 I  BIT  #200,@#I0C  ; TEST  DONE  BIT 

BEQ  L00P2 

MOV  R 1,0# IOC  :S£T  ADDRESS  TO  CH  NCI 
MOV  0#IOI,(R4)+  J  READ  CH  NC2 
DEC  RO 

BNE  TEST 1  ; CHECK  to  see  IF  N  SAMPLES  HAVE  BEEN  TAKEN 

RTS  PC  JIF  SO,  RETURN 

.END 


Enclosure  2  to  Appendix  B 


n  isv  e  i  oc^Mcn»  >a  *<p<  •  on 

r 

•'  •jiioopmt »►(£  cno  'zQL,-> I K,,3  5 1 rA'llEOLlS  L  I NIEA*  ALGE?°(I"'IC 

C  EGUATIONS  USING  GAUSS  ELINIMATION  '4lTH  PGM  "INGOING. 

C 

C  eon*<  "CIPC'J I T  T1-ISCPY.  A  C0NP,2TATI0NAL  APoo0ACH“  SY  S.  14. 

C  0 IPECTGP  t  PAGE  211. 

C 

C  FOR*  OF  CALL  IS! 

C  CALL  GAUSS 

C  THE  NYM  NATRIY  A  5r,JE  M  "EC‘rOP  B ,  THE  N  '.'SCTOP  Y  AMS  TUE  D I YENS I ON 

C  N  HtJST  SE  COHHGN  UAPIABLES. 

c 

subroutine  gauss 

CONNOM  A  <44, 44),  3  (44)..  m 

DO  5  I -1  ,M 
1 1  -*  I  ♦  1 

IF  ( ABS ( A<  I  ••  I  >  >  •  L£ .  1 . 5— 10>  GO  T0  1 
GO  TO  if 

1  CONTINUE 

IF  (I.EG.N)  GO  T0  10 
DO  I*1  J  •»  1 1 .  N 

IF  > ABS< A< J . I > > .LS. : .E-10?  GO  T0  14 
I  0 1 J 
GO  T0  IS 
14  CONTINUE 

GO  TO  10 

16  DO  2  K-l.M  . 

PIV*A< IPIUrK >  »  * 

_  A  (  IPI".  K  )  "  A  <  I  .  r.  ) 

2  A  <  I  ,  K  ?  ■*  p  I 
pi'.'-e<  1  p  1  '.*  > 

3 <  IPI'.'}  »e<  i  > 

3  t  r  \  -o  tii 

If  IF  (I.EQ.N?  GO  T0  2 

DO  9  J  1  *  1 1 . N 

8  rV  I  J1  ?  -A<  I  J1  >  AV  I  .  I  > 

9  '  I  >  "9  •  I  '  A  *  I .  I } 

DO  f  J  “  1 1 .  N 
DO  4  K -  1 1 . M 

4  A  *  J  K  '  *  A  *  J  •  '  —A  1  j  ,  I )  —A  >  I .  v  ) 

f  3  <  J  '  "3  <  J  '  -9  ‘  I  '  -A  <  .*  I  / 

2  2  «.  N  >  ’9  • N •  A  <  N  .  N  :• 

r.n  c  V~2  ■*’ 


DO  ?  J  ->L  .  N 

?  SUM-SUNfA*:  I  .  J>~9<  J? 

G  3  < I } *9 • I ? -SUN 

GO  tq  n 

10  UP  I TE ? .  2  > 

9  FORNAT (  '  EQUATIONS  ARE  LINEARLY  DEPENDENT'? 

STOP 

1 1  RETURN 
END 


u  u  u 


1325 

1950 

2000 

2100 

2200 

2300 

3000 
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DO  1925  ILQOP*  1  ,NLQ0P 
CONTINUE 

PLOT  OUT  DATA 

CALL  DAOUT  ( IX 1 ( I+N/2 > , INT < OUT ) > 

FORMAT  ( 217 »  F10 . 2 ) 

CONTINUE 

WRITE  (7,2100)  AUG/ (NS-NT+1 ), AIN/AUG 

FORMAT  ('  AVERAGE  OUTPUT  SQUARED ' ,F14.2 , '  SNR  IMPROVEMENT' 
WRITE  (7,2200) 

FORMAT  ('  NEW  DATA,  SAME  WEIGHTS?  1=YES,  0-  NO') 

READ  (5,2300)  I FLAG 
FORMAT  (16) 

IF  ( I FLAG. EG. 0)  GO  TO  1050 
CALL  XTSMP2 ( N , NC 1 ,NC2» I X 1 , 1X2) 

GO  TO  1600 

STOP 

END 
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1000  CONTINUE 

1050  WRITE  (7,11 00 ) 

1100  FORMAT  ('  NUMBER  OF  WEIGHTS , NUMBER  OF  DELAY  LOOPS') 

READ  (5.1200)  N.NLOOP 
1200  FORMAT i  216 ) 

IF  (N.GT.70)  GO  TO  3000 

C 

C  CALCULATE  A  it  B  FOR  GAUSS  SUBROUTINE 

C 

DO  1300  J  *  1  ,  N 

DO  1225  K  *  1 »  N 

IF  ( IFLAG . EG . 0 )  GO  TO  1240 
ID= IASS ( J-« ) f 1 
A< J,K)*C( ID) 

1225  WRITE  (7,1230)  J,K,A(J,K) 

1235  FORMAT  ('  A(  ' , 12,  ' ,  '  . 12,  ' )  *  ',E12.5) 

1240  B(J)-0.0 

DO  1250  1*1,  NS-*’  1  -NT 

1250  8 ( J ) *8 ( J ) ♦FLOAT ( 1X2 ( If J-l ) >*FLOAT( IX 1 ( IfN/Z) ) 

IF  (IFLAG. EG. 0)  GO  TO  1300 
WRITE  (7,1280)  J , 8 ( J ) 

1280  FORMAT  i'  B(',I2,')  *  ',E12.5) 

1300  CONTINUE 

C 

C  USE  GAUSS  ELIMINATION  SUBROUTINE  TO  CALCULATE  OPTIMUM  WEIGHTS 

C 

CALL  GAUSS 
C 

C  PRINT  OUT  WEIGHTS 

C  * 

DO  1400  J  *  1 , N 

WRITE  (7,1500)  J  , 8 ( J  ) 

1400  CONTINUE 

1500  FORMAT  iIG,F13.7) 

1600  AUG *0.0 

C 

C  CALCULATE  FILTER  OUTPUT  USING  OPTIMUM  WEIGH~S 

C 

DO  2000  I*1,NS+1-NT 

OUT “FLOAT ( 1X1 ( IfN/Z) ) 

DO  1800  J  *  1 , N 

OUT *OUT-B ( J ) *  1X2 ( I  +  J-l  ) 

1900  CONTINUE 

C 

C  CALCULATE  OUTPUT  POWER 

C 

AVG*AUGfOUT**2 

C 

C  PRINT  OUT  INPUT  AND  OUTPUT  DATA 

C 

WRITE  (7,1350)  I , 1X1 ( If N/Z > , OUT 
C 

C  DELAY  SO  THAT  STRIP  CHART  RECORDER  CAN  KEEP  UP 

C 
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PROGRAM  OPTFIL 

C  B.  8.  PETERSON  10  NOV  S3 

C  PROGRAM  READS  IN  TWO  CHANNELS  OF  ANALOG  DATA  USING  SUBROUTINE 

C  XTSMP2 .  THE  OPTIMUM  WEIGHTS  OF  FIR  FILTERS  OF  VARIOUS  LENGTHS 

C  ARE  CALCULATED  USING  A  ONE  SHOT  LINEAR  LEAST  SQUARES  TECHNIQUE. 

C  THE  AVERAGE  POWER  IN  THE  OUTPUT  IS  THEN  CALCULATED  FOR  EACH  LENGTH. 

C 

c  arrays: 

C  1X1  &  1X2 

C  A  &  B 

C 
C 
C 

C  C 

C 
C 

INTEGER  IXl <3000) , 1X2(3000) 

COMMON  A ( 30 , 30  > , 8 ( 30 ) , N 
REAL  C ( 30 ) 

C 

C  INPUT  PARAMETERS 

C 

30  WRITE  (7,100) 

100  FORMAT  ('  ENTER  #  SAMPLES , S+N  CH , NOISE  CH, LARGEST  F I LTER  LENGTH') 

READ  (5,400)  NS , NC 1 , NC2 , NT 
400  FORMAT  (417) 

WRITE  (7,300) 

500  FORMAT  ('  type  DAT  A"7  ( 1 = YES , 0  =  N0 >  '  ) 

READ  ( 5 , SOO  )  I FLAG 
BOO.  FORMAT  (14) 

A  I N  =  0 . 0 
DO  700  J  3 1 , NT 
G30  C(J)=0.0 

700  CONTINUE 

C 

C  DELAY  BEFORE  SAMPLING 

C 

DO  S50  1=1 , 10000 
S30  CONTINUE 

N1 =23G*NC1+16 
N2=23G*NC2>16 

CALL  XTSMP2(NS,N1 ,N2, 1X1 , 1X2) 

IF  (IFLAG.EQ.O)  GO  TO  940 
900  DO  920  1=1 ,NS 

920  WRITE  (7,930)  I , 1X1 ( I ) , 1X2 ( I ) 

930  FORMAT  (317) 

940  DO  1000  1=1, NS+ 1 —NT 

930  DO  990  J  =  1  ,  NT 

C( J ) *C( J )+FLOAT( IX2( I+J-i ) ) *FLOAT ( 1X2 ( I ) > 

990  CONTINUE 

C 

C  CALCULATE  INPUT  POWER 

C 

AIN=AIN+FLOAT ( 1X1 ( I ) ) *FLOAT ( IXl ( I ) ) 


INPUT  DATA  FROM  TWO  MICROPHONES 
MATRIX  AND  VECTOR  IN  SOLUTION  OF 

SpASS^flClieVSREg^T«gS^^RsiD92NSuISe 

IN  THE  SOLUTION  OF  AX=9) 

AUTOCORRELATION  VECTOR  OF  THE  REFERENCE 
INPUT  (1X2)  USED  TO  CALCULATE  A. 


However, in  the  program  (enclosure  1)  because  of  the  long  time 
required  for  the  large  number  of  calculations,  they  were  set 
equal  which  allowed  the  calculation  of  only  one  autocorrelation 
vector  (C).  The  longer  the  data  set  the  more  valid  this 
approximation  is. 

After  the  A  matrix  and  the  B  vector  are  calculated,  the 
optimum  weight  vector  is  calculated  using  a  Gauss  elimination 
subroutine  (enclosure  2).  The  filter  output  was  then  calculated 
using  this  weight  vector  and  the  output  energy  compared  to  the 
primary  input  energy. 

The  results  were  not  impressive.  Data  from  both  the  Luder 
yawl  engine  and  the  gasoline  research  were  analyzed  with  minimal 
reduction  in  noise  power.  In  view  of  these  results  and  others 
obtained  since,  the  explanation  of  the  poor  performance  of  the 
adaptive  early  filters,  was  there  limited  length  and  not 
necessarily  slow  convergence. 


Repeated  versions  of  equation  B1  can  be  written  in  matrix  form: 


Y  =  P  -  RH  ( B2 ) 

where  POO  =  x^k-N/2)  and  Rki  =  x2(k-i).  The  total  output 
energy  is  given  by: 


E  = 


YTY 


=  PTP  - 


HTRTP 


htrtrh 


(B3) 


This  energy  is  minimized  by  setting  it's  gradient  with  respect  to 
H  equal  to  0. 

VhE  =  -2  RTP  ♦  2  RTRH  =  0  ( B4 ) 

or  AH  =  B 

where  A  s  RTR  and  B  =  R^P. 

The  elements  of  A  and  B  are  given  by: 

Aij  =  Aji  S  XgOc-i)  x2(k-j) 

Data  Set  (B5) 

Bj  s  x2(k+j)  x1 ( k+N/2 ) 

Data  Set 

Essentially  A  is  an  N  x  N  autocor relation  matrix  of  the  reference 
input  and  B  is  a  M  vector  of  the  cross-correllation  of  the 
reference  input  and  the  delayed  primary  input.  Strictly  speaking 
A .j  i  ft  A22  and  A12  ft  A^  etc.  because  A^  contains  one  early 


sample  not  in  A,,  and  A??  contains  a  late  sample  not  in  A... 
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FORTRAN  Program  to  Calculate  Optimum  Filter  Weights 

In  the  LMS  adaptive  noise  cancellation  scheme  the  weights  of 
the  adaptive  filter  are  adjusted  to  minimize  the  output  power. 

In  previous  FORTRAN  implementations  of  the  LMS  algorithm,  the 
filter  was  found  to  ineffective  in  canceling  engine  noise.  It 
was  felt  this  may  be  due  to  the  weights  not  converging  within  the 
space  of  the  data  set.  Using  off-line  parameter  estimation 
techniques  it  is  possible  to  calculate  the  optimum  fixed  weight 
filter  for  a  given  data  set,  and  then  analyze  this  filter  using 
the  data  set.  The  filter  is  optimum  in  the  sense  that  the 
remaining  output  power  for  the  data  set  is  the  minimum  of  all 
possible  weight  vectors.  The  filter  is  illustrated  in  Figure  B1. 


For  a  fixed  weight  vector  H,  the  filter  output,  y(k)  is  given 


by : 

y(k) 


N- 1 

x^k-N/2)  -  ^  hi  XgOc-i) 
UO 


CB 1 ) 


B1 


The  monitor  on  the  TMS-32Q  development  system  allows 
assembly  language  programs  to  be  down-loaded  from  another 


C5 


ti 


I 

i 


r 

i 


computer.  It  then  assembles  the  program  into  its  machine  code  for 
execution.  To  analyze  the  effects  of  various  numbers  of  weights, 
sampling  frequencies,  and  adaptive  gain  coefficients  would 
require  many  modifications  to  the  assembly  language  program.  To 
make  these  modifications  easy,  a  BASIC  7  program  was  developed  on 
the  Dartmouth  Time  Sharing  system  to  generate  TMS-320  assembly 
code.  The  program  is  interactive  in  that  it  asks  for  the  number 
of  weights,  sampling  frequency,  and  adaptive  gain  coef f icie  it.  It 
produces  straight  line  code  so  the  time  spent  on  the  tapped  delay 
line  filtering  and  weight  adaptation  is  minimized.  There  are 
delay  loops  added  to  the  code  to  provide  a  range  of  sampling 
frequencies.  Enclosure  Cl  is  a  listing  of  the  Basic  program. 
Enclosure  ’C2  is  a  listing  of  the  assembly  language  code  produced 
by  this  program  for  a  10  weight  filter. 

This  system  has  been  implemented  to  cancel  engine  background 
noise  in  voice  communications.  Table  Cl  summarizes  the  amount  of 
noise  reduction  for  various  numbers  of  weights,  sampling 
frequencies,  and  adaptive  gain  coef f icients .  The  amount  of  noise 
reduction  is  calculated  by  comparing  the  RMS  voltage  of  the 
3ignal  plus  noise  input  to  the  RMS  voltage  of  the  output.  Anti¬ 
aliasing  filters  set  at  half  the  sampling  frequency  were  used  on 
both  input  signals. 


C6 


From  the  data  shown  in  Table  Cl ,  it  is  clear  that  the  amount 
of  noise  reduction  is  improved  by  increasing  the  number  of  filter 
weights.  At  the  higher  sampling  rates,  there  in  very  little  noise 
reduction.  At  the  lower  sampling  frequencies  the  noise  reduction 
is  even  more  apparent  (  8  dB  at  2  kHz.).  This  indicates  that  to 
achieve  the  same  amount  of  noise  reduction  at  a  sampling 
frequency  of  10  kHz  would  require  five  times  more  weights  or  on 
the  order  of  300  weights.  Also  there  is  an  apparent  increase  of 
signal  reduction  when  the  adaptive  gain  coefficient  is  increased. 
This  is  not  only  an  increase  in  noise  reduction;  but,  because  the 
filter  is  adapting  so  quickly  to  the  signal,  it  is  also  canceling 
some  of  the  desired  signal  and  not  just  the  noise. 

The  noise  cancellation  system  only  performed  marginally  in 

% 

the  engine  noise  application  at  useful  sampling  frequencies.  To 
improve  the  performance  would  require  a  processor  capable  of 
implementing  more  weights.  The  next  generation  of  TMS-320  may 
have  this  capability.  This  work  has  shown  that  adaptive  filtering 
with  a  maximum  of  68  weights  is  realizable  for  speech  processing 
and  may  have  other  applications. 


Sampling  Number  of 
Frequency  KHz  Weights 


Adaptive 

Gain 


Signal 

Reduction  d3 


10 

68 

14 

3.5 

10 

68 

12 

2.6 

10 

68 

10 

2.4 

10 

68 

8 

1.4 

10 

68 

6 

* 

10 

50 

14 

3.2 

10 

50 

12 

2.3 

10 

50 

10 

2.3 

10 

50 

8 

1.5 

10 

50 

6 

« 

10 

30 

14 

2.5 

10 

30 

12 

1.9 

10 

30 

10 

1.9 

10 

30 

3 

1.5 

10 

30 

6 

• 

5 

68 

14 

2.4 

5 

68 

12 

2.4 

5 

68 

10 

2.1 

5 

68 

8 

-.8 

5 

68 

6 

* 

5 

50 

14 

2.1 

5 

50 

12 

2.2 

5 

50 

10 

2.2 

5 

50 

8 

0.0 

5 

50  ‘ 

6 

• 

5 

30 

14 

1.9 

5 

30 

12 

1.9 

5 

30 

10 

1.9 

5 

30 

8 

0.9 

5 

30 

6 

• 

2 

68 

14 

8.4 

2 

68 

12  . 

7.9 

2 

68 

10 

6.0 

2 

68 

8 

0.6 

2 

68 

6 

• 

2 

50 

14 

6.4 

2 

50 

12 

5.6 

2 

50 

10 

4..  9 

2 

50 

8 

0.4 

2 

50 

6 

-.4 

2 

30 

14 

2.4 

2 

30 

12 

2.4 

2 

30 

10 

2.1 

2 

30 

8 

-.8 

2 

30 

6 

» 

(*  -  Did  not 


track) 


Table  Cl.  Test  results  of  Noise  Cancellation  System 
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100  ! 

110  !  PROGRAM  NAME:'  ADAPT-F  (ADAPTIVE  FILTER) 

120  !  VERSION  4 

130  ! 

140  !  PROGRAMMER:  K.  U.  DYKSTRA 
150  ! 

160  !  DATE:  20  SEPT  84 
170  ! 

180  !  THIS  PROGRAM  PRODUCES  THE  ASSEMBLY  LANGUAGE  ADAPTIVE 
190  !  FILTER  PROGRAM  FOR  THE  TMS-320  EVALUATION  MODULE.  IT 
200  !  CAN  BE  USED  IN  CONJUNCTION  WITH  THIS  MODULE  TO  DOWN 

210  !  LOAD  THE  ASSEMBLY  LANGUAGE  PROGRAM.  A  TAPPED  DELAY 

220  !  LINE  ALGOR I TM  IS  USED  FOR  THE  ACTUAL  FILTERING.  THE 
230  !  WEIGHTS  ARE  ADJUSTED  USING  A  LMS  ALGORITHM.  INPUTS 
240  !  ARE  NUMBER  OF  WEIGHTS,  SAMPLING  FREQ,  ADAPTIVE 
COEFFICIENT. 

250  !  THE  CODE  THAT  IS  PRODUCED  IS  STRAIGHT  LINE  CODE  WITH 
260  !  NO  LOOPING  FOR  MAXIMUM  SPEED. 

270  ! 

230  !  DATA  MEMORY  ASSIGNMENT 

290  ! 

300  1  0-68  REFERENCE  INPUTS 

310  69  CHANNEL  1  CONTROL 

320  !  70  CHANNEL  2  CONTROL 

330  !  71  MASK  FOR  INPUTS 

340  !  72  MASK  FOR  OUTPUTS 

350  !  '73  FILTER  OUTPUT 

360  !  74  SIGNAL  +  REFERENCE  INPUT 

370  !  75  OUTPUT 

380  !  76-143  FILTER  WEIGHTS 

390  ! 

400  !  INPUT  NUMBER  OF  WEIGHTS. 

410  ! 

412  DO 

4  -|  4  let  ERROR  s  0 

420  PRINT  "ENTER  NUMBER  OF  WEIGHTS  (1-68)." 

430  INPUT  NWTS 

432  IF  NWTS  <  1  OR  NWTS  >  68  THEN  LET  ERROR  =  1 
434  LOOP  UNTIL  ERROR  =  0 
440  ! 

450  !  INPUT  SAMPLING  FREQUENCY  . 

460  ! 

470  LET  LOWF  =  1  /  ( .  00Q2*(37+( 6 *NWTS )+(20*256)  ) ) 

472  LET  HIGHF=1/( . 0002* ( 37+ ( 6 »NWTS ) +20 ) ) 

474  IF  HIGHF  >  15.6  THEN  LET  HIGHF  =  15.6 
476  DO 

473  LET  ERROR  =  0 

480  PRINT  "ENTER  SAMPLING  FREQUENCY  (";LOWF;  "  TO  "; HIGHF;" 
KHZ )  " 

490  INPUT  FREQ 

495  IF  FREQ  <  LOWF  OR  FREQ  >  HIGHF  THEN  LET  ERROR  =  1 

496  LOOP  UNTIL  ERROR  =  0 


500  LET  DELAY = I NT ( C  C 1 /(FREQ*. 0002) )-33-(6*NWTS) )/2Q) 

510  LET  AFREQ=1/C .0002*(37+(6*NWTS)+(20*DELAY) ) )  ' 

520  PRINT  "SAMPLING  FREQUENCY  =  " ; AFREQ; ”  DELAY  =  " ; DELAY 
530  ! 

5  32  DO 

534  LET  ERROR  -  0 

540  PRINT  "ADAPTIVE  COEFFICIENT  s  1  / ( 2** ( 1 6-N ) ) ,  ENTER  N” 
550  PRINT  "IF  YOU  WANT  TO  DOWN  LOAD  END  WITH  <CTRL  C>  ELSE 
<RETURN>" 

560  INPUT  COEF 

562  IF  COEF  <  0  OR  COEF  >  15  THEN  LET  ERROR  =  1 
564  LOOP  UNTIL  ERROR  =  0 
570  ! 

580  !  SET  UP  DATA  CONSTANTS 
590  ! 


600 

PRINT 

It  >  ft 

610 

PRINT 

If 

AORG  0" 

620 

PRINT 

If 

B  BEGIN" 

630 

PRINT 

If 

NOP" 

640 

PRINT 

ft 

NOP" 

650 

PRINT 

If 

DATA  >7FFQ" 

660 

PRINT 

It 

DATA  >8000" 

670 

PRINT 

"BEGIN 

LACK  4" 

630 

PRINT 

If 

TBLR  71" 

690 

PRINT 

If 

LACK  5" 

700 

PRINT 

It 

TBLR  72" 

710 

PRINT 

It 

LACK  3" 

720 

PRINT 

It 

SACL  69" 

730 

PRINT 

ft 

LACK  >83" 

740 

PRINT 

ft 

SACL  70" 

750 

PRINT 

ft 

SOVM" 

760 

PRINT 

It 

LARP  1" 

770 

780 

j 

!  INPUT  REFERENCE  INPUT 

790 

800 

i 

PRINT 

"LOOP 

OUT  69,0" 

810 

PRINT 

It 

NOP" 

820 

PRINT 

If 

NOP" 

830 

PRINT 

If 

IN  0,2" 

840 

PRINT 

ft 

ZALS  71" 

850 

PRINT 

It 

XOR  0" 

360 

PRINT 

It 

SACL  0" 

370 

380 

J 

!  DELAY  FOR  SAMPLING  FREQUENCY 

890 

900 

i 

PRINT 

If 

LARP  0" 

910 

PRINT 

It 

LARK  0,"; DELAY 

920 

PRINT 

ffD1 

NOP" 

921 

FOR  I 

=  1  TO 

7 

922 

PRINT 

it 

NOP" 

923 

NEXT 

L 

930 

PRINT 

ft 

BANZ  D 1 " 

940  ! 

950  !  PERFORM  TAPPED  DELAY  LINE  FILTER 


960  ! 

970  PRINT  "  ZAC’* 

975  PRINT  '•  MPYK  O'* 

980  PRINT  ’’  LARP  1” 

990  PRINT  "  LARK  1,143" 

1000  FOR  I  =  NWTS  -  1  TO  0  STEP  -1 
1010  PRINT  "  LTD  " ; I 

1020  PRINT  "  MPY 

1030  NEXT  I 

1040  PRINT  "  APAC" 

1050  PRINT  "  SACH  73” 

1060  ! 

1070  !  OUTPUT  =  SIGNALAREF  -  FILTERED  OUTPUT 
1080  ! 

1090  PRINT  "  LAC  74" 

1100  PRINT  "  SUB  73” 

1110  PRINT  "  SACL  75” 

1120  ! 

1130  !  COEF*OUTPUT  >  T 
1140  ! 

1150  PRINT  "  LAC  75  ,  '* ;  COEF 

1160  PRINT  "  SACH  73” 

1170  PRINT  "  LT  73" 

1180  ! 

1190  !  OUTPUT 

1200  ! 

1210  PRINT  "  ZALS  72"  . 

1220  PRINT  "  XQR  75" 

1230  PRINT  "  SACL  75" 

1240  PRINT  "  OUT  75,2" 

1250  ! 

1260  !  INPUT  REFERENCEASIGNAL 

1270  ! 

1280  PRINT  "  OUT  70,0" 

1290  PRINT  "  NOP" 

1300  PRINT  ”  NOP" 

1310  PRINT  ”  IN  74,2" 

1320  PRINT  "  ZALS  71" 

1330  PRINT  "  XOR  74" 

1340  PRINT  "  SACL  74" 

1350  ! 

1360  !  DELAY  FOR  SAMPLING  FREQUENCY 

1370  ! 

1380  PRINT  "  LARP  0" 

1390  PRINT  "  LARK  0,"; DELAY 

1400  PRINT  "D2  NOP" 

1401  FOR  I  =  1  TO  7 

1402  PRINT  "  NOP" 

1403  NEXT  I 

1410  PRINT  "  BANZ  D2” 

1420  PRINT  "  LARP  1" 

1430  PRINT  "  LARK  1 , 143" 

1440  ! 

1450  !  PERFORM  ADAPTIVE  PART 


1460  ! 

1470  FOR  I  =  NWTS-1  TO  0  STEP  -1 
1480  PRINT  "  ZALH  *" 

1490  PRINT  "  MPY  " ;  1+1 

1500  PRINT  "  APAC" 

1510  PRINT  "  SACH 

1520  NEXT  I 

1530  PRINT  "  8  LOOP" 

1540  PRINT  "  END" 

1550  PRINT  "<” 

1560  END 
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ENTER  NUMBER  OF  WEIGHTS  (1-68). 

?  10 

ENTER  SAMPLING  FREQUENCY  (  .958405  TO  15.6  KHZ). 

?  10.0 

SAMPLING  FREQUENCY  =  10.0604  DELAY  =  20 

ADAPTIVE  COEFFICIENT  =  1 /(2**( 1 6-N) ) ,  ENTER  N 
IF  YOU  WANT  TO  DOWN  LOAD  END  WITH  <CTRL  C>  ELSE  <RETURN> 
?  12 
> 

AORG  0 
B  BEGIN 
NOP 
NOP 

DATA  >7FFO 
DATA  >8000 
BEGIN  LACK  4 
‘ TBLR  71 
LACK  5 
TBLR  72 
LACK  3 
SACL  69 
LACK  >83 
SACL  70 
SOVM 
LARP  1 

LOOP  OUT-  69,0 
NOP 
NOP 

IN  0,2 
ZALS  71 
XOR  0 
SACL  0 
LARP  0 
LARK  0,  20 
D 1  NOP 
NOP 
NOP 
NOP 
NOP 
NOP 
NOP 
NOP 

BANZ  D 1 
ZAC 

MPYK  0 
LARP  1 
LARK  1,143 
LTD  9 
MPY  »- 
LTD  3 
MPY  »- 
LTD  7 


MPY  »- 
LTD  6 
MPY  »- 
LTD  5 
MPY  *- 
LTD  4 
MPY  •- 
LTD  3 
MPY  »- 
LTD  2 
MPY  •- 
LTD  1 
MPY  *- 
LTD  0 
MPY  *- 
APAC 
SACH  73 
LAC  74 
SUB  73 
SACL  75 
LAC  75,  12 
SACH  73 
LT  73 
ZALS  72 
XOR  75 
SACL  75 
OUT  75,2  •  • 
OUT  70,0 
HOP 
MOP 

IN  74,2 
ZALS  71 
XOR  74 
SACL  74 
LARP  0 
LARK  0,  20 
D2  NOP 
NOP 
NOP 
NOP 
NOP 
NOP 
NOP 
NOP 

BANZ  D2 
LARP  1 
LARK  1,143 
ZALH  » 

MPY  10 
APAC 
SACH  *- 
ZALH  * 

MPY  9 
APAC 


SACH  »- 
ZALH  * 
MPY  3 
APAC 
SACH  »- 
ZALH  * 
MPY  7 
APAC 
SACH  »- 
ZALH  » 
MPY  6 
APAC 
SACH  *- 
ZALH  * 
MPY  5 
APAC 
SACH  *- 
ZALH  * 
MPY  4 
APAC 
SACH  «- 
ZALH  * 
MPY  3 
APAC 
SACH  *- 
ZALH  * 
MPY  2 
APAC 
SACH  »- 
ZALH  » 
MPY  1 
APAC 
SACH  »- 
B  LOOP 
END 
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Notch  Filter  Analysis  and  Results 

In  [1]  and  [2]  it  was  pointed  out  that  when  the  reference 
input  (x2(k)  in  Figure  D1)  is  periodic  in  the  length  N,  of  the 
adaptive  filter,  the  result  is  a  notch  filter  from  x-,00  to  y(k) 


Figure  D1 

Because  the  background  noise  due  to  an  engine  is  quite  periodic, 
it  follows  that  any  reference  that  contains  all  the  harmonics  of 
the  engine  noise  should  suffice.  A  periodic  reference  input 
consisting  of  one  "I"  followed  by  N-1  "Q"'s  i.e. 

x2(k)  s  1  for  k  =  mN  (D1 

=  0  for  k  =  mN 

where  M  is  both  the  number  of  samples  per  period  of  the  engine 
noise  and  the  number  of  filter  weights  and  m  is  an  integer  not 


D2 


only  contains  all  these  harmonics  but  also  makes  possible 
implementing  a  filter  of  large  N  in  real  time. 


The  filter  and  adaptive  algorithm  equations  now  become: 

N-1 

y(m N+i)  =  x^mN+i)  -  h^(mN>i)  x2(mN+i-j) 

1=0 

=  x1 (mN+i)  -  h^mN+i) 


(D2) 


and 


ht(  (m+1  )N-fi)  =  h.(mN+t)  +  a  y(k)  (D3) 

Because  the  algorithm  requires  only  one  subtraction,  one  addition 
and  one  multiplication  per  sample  period,  it  can  be  implemented 
very  efficiently  with  the  order,  N ,  limited  only  by  memory. 
Further,  if  a  s  2“n  the  multiply  can  be  accomplished  with  right 
shifts  allowing  implementation  on  a  general  purpose 
microprocessor . 

Equations  (D2)  and  C  D  3 )  can  be  written  using  z  transforms: 

Y(z)  =  X1 ( z )  -  H.(z)  ( D4 ) 

zN  H^z)  =  Hi(z)  -  a  Y(z)  (D5) 

Combining  (D4)  and  (05)  results  in  the  transfer  function  from 
x 1 ( k)  to  y( k) : 

Hz)  =  -  1 

X1  ( z)  zN  -  ( 1 -a)  (06) 

The  N  zeros  of  (06)  are  on  the  unit  circle  at  eJ^",r'</N,  k  =  0  to  N- 
1  and  each  zero  has  a  corresponding  pole  at  (1-a)1/N  A 

pole-zero  plot  for  a  =  0.5  and  N  =  20  is  shown  in  Figure  D2.  The 
actual  filter  implemented  had  a  =  0.25  and  N  =  400.  Each  pole 


was  therefore  at  a  radius  of  0.999231.  Figure  D3  is  a  small 
section  of  the  pole-zero  plot  for  this  case  illustrating  the 
seven  poles  and  zeros  nearest  z=1.  The  net  effect  is  a  notch 
filter  with  a  notch  at  each  harmonic  of  x^OO.  The  width 
(distance  between  -3  db  points)  of  each  notch  is  approximately: 

a  f  a  f 

s  _  o 

where  f  and  f  are  the  sampling  and  fundamental  frequencies 
respectively.  For  N  =  400,  as  0.25  and  f  s  7  kHz,  the  notches 
are  1.4  Hz  wide  and  adjacent  notches  are  separated  by  17.5  Hz. 
This  frequency  response  is  impossible  to  verify  using  our  analog 
spectrum  analyzers.  Figure  D4  shows  the  output  spectrum  for  N  = 
160,  a  =  0.25  and  f  s  40  kHz.  The  input  was  a  sinusoid  swept 
from  20  Hz  to  10  kHz.  The  notches  are  19.9  Hz  wide  and  spaced  at 
250  Hz  intervals.  The  slight  rolloff  at  high  frequencies  is  due- 
to  the  hold  function  of  the  D/A  converter. 

If  the  output  of  the  filter  is  taken  at  the  output  of  the 
adaptive  filter  before  the  summer,  the  resulting  filter  is  a 
recursive  comb  filter  passing  only  those  components  of  x^Ck)  that 
are  periodic  in  N  samples.  Similar  non-r ecur s ive  comb  filters 
that  exploited  the  periodicity  of  the  speech  waveform  to 
eliminate  white  background  noise  were  proposed  in  [3]  and  [4], 
These  filters  were  adaptive  in  the  sense  that  the  period  of  the 
filter  was  continuously  adapted  to  match  the  period  of  the  speech 
waveform.  This  type  of  filter  is  also  the  digital  version  of  the 
analog  waveform  eductors  of  many  years  ago. 

If  x  £  ( >< )  is  an  arbitrary  periodic  signal  with  a  period  of  M 
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°RCGRAM  ADCDLS 
3.  3.  PETSRSON 

Q°OGRAM  TO  IMPLEMENT  ADAPttmE  CI1TEP  IN!  Ft5EGLlcl'lCY  nOMA*N 
USING  WEIGHTED  1£AST  SGUARES  ADAPTIVE  ALGORITUM 
REF1  M.  DENTIMO,  J.  MCCOOL,  1  3.  WIDROW,  "ADAPTIVE  "ILTERING 
IN  THE  CREQUENCY  DOMAIN",  PRCC  IEEE/  'CL.  EE /NO.  12.  DEC  "3. 
MAXI NUN  SIZE  OF  EACH  BLOCK  OF  data  IS  712 

D I  MENS  I  ON  X 1  ( 227  >  ,  Y 1  (  27T  >  ,  X2  (  227  )  Y2  ( 227  )  .•  YO  (  227  )  .•  YO  (  227 ) 
DIMENSION  WR ( 712) /WI (712? rXH(27S? , YH<276) »  P2 ( 22S ? , CR(27S) 
DIMENSION  Cl (226) 

INTEGER  II  (7120) , 12(7120) 


VARIABLES 

XI  i  Y 1  PRIMARY  INPUT  (S0-rH  TIME  AND  FREOUEMCY  DOMAIN  DATA) 
;<2  1  Y2  REFERENCE  INPUT 
XO  &  YO  OUTO,JT 

YH  £  YH  WEIGHTS  OF  ADAPTIVE  FILTER  < FREQUENCY  DOMAIN? 

WR  £  WI  TABLES  OF  COSINES  AND  SINES 

11  ttmc  hqmain  po//«APY  INPUT 

12  T I ME  DOMAIN  occcspijc  rn|piiT 

CALCULATE  vectors  of  sines  and  cosines 

NMAX*7 1 2  '  •  * 

TPQN*S . 2321S/eL0AT ( NMAX)  *  * 

DO  3  I*l>  NMAX 

OU-CI_0AT  (  <  l-l  )  )  *T  3Qh| 

WR ( I ) *COS( ou ) 

W I t  I ) 3-S I N ( ou  » 

continue 
input  °apameters 

WRITE  '7,10' 

FORMAT  <  '  NEW  DAT(l=Y),2ER  Wt,#  SAMP. 8LK  SZ.CH  #S,ADP  GAIN') 
READ  '7/20?  IF1 , IF2/N,MW,NC1 /NC2.G 
FORMAT  (SI6/F10.8) 

NW=2**tw 
OMG  *  1 . 0-G 

IF  (IF1.ME.1)  GO  T0  100 
N  l  =  276 *NC  1  *  IE 
N2  *  27SJ»NC2  +  IS 

INSERT  DELAY  9E-0RE  SAMPLING 

DO  30  J  = 1 , 20000 

continue 

CALL  X T S M P 2  (  N  ,  Ml  ,  N2 , 1 1  ..  12  ) 

NW2  *  NW/2 

IF  < IF2. EG. 0?  GO  tq 


results  are  subtracted  from  X^f),  giving  the  FFT  of  the  output, 
which  is  then  inverse  transformed  to  realize  the  time  domain 
output.  The  data  is  processed  in  blocks  of  N  (=2n)  points. 

Since  each  complex  element  of  Y(f)  is  a  function  only  of  the 
corresponding  coefficient  of  H(f)  and  not  the  entire  vector  the 
adaptive  algorithm  can  be  made  to  have  much  more  efficient  and 
predictable  convergence  properties.  Therefore,  unlike  the  time 
domain  problem,  when  doing  off-line  processing  on  the  LSI-11/2 
with  long  filter  lengths  but  limited  data  storage,  in  the 
frequency  domain  filter  one  can  be  assured  the  weights  will 
converge  within  the  data  set.  Also,  because  of  the  efficiency  of 
the  algorithm,  the  data  can  be  processed  in  reasonable  time. 

In  the  FORTRAN  implementation  (enclosure  1)  the  complex  LMS 
algorithm  proposed  in  [6]  has  been  changed  to  a  weighted  least 
squares  algorithm  -for  more  efficient  and  predictable  convergence 
properties.  Each  element  of  H(f)  is  the  best  least  squares 
estimate  with  the  input  data  exponentially  weighted.  The  filter 
can  therefore  track  time  varying  parameters  with  the  time 
constant  of  the  exponential  weighting. 
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Frequency  Domain  Adaptive  Filtering 

A  fundamental  problem  in  both  real  time  implementation  and  in 
off-line  analysis  of  adaptive  noise  cancellation  is  that  the 
filter  length  is  severely  limited  by  memory,  processor  speed,  or 
both.  At  the  expense  of  much  more  complicated  software,  digital 
filters,  including  adaptive  filters,  can  be  implemented  more 
efficiently  in  the  frequency  domain  than  the  time  domain.  The 
basic  structure  of  a  frequency  domain  adaptive  noise  canceler, 
proposed  in  [5]  is  shown  in  Figure  El. 


Figure  El 

The  Fast  Fourier  Transforms  (FFT's)  of  the  two  signal  inputs 
are  calculated,  the  complex  coefficients  of  the  reference  input 
FFT  are  multiplied  by  complex  filter  coefficients,  H(f).  The 


E. 
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.title  '"alt: 

. GLOSL  REAuTI 

3.  3.  PETERSON,  3  MAY  34 

IMPLEMENT’S  NOTCH  FIL”ER  IN  REAL  ”IME. 

CALLED  FROM  FGF”RAN  BY  "CALL  REALTI  (NC.nw) 

NC  *  253+  CH .  NC .  +  xE 

NU  =  NO.  WEIGHTS  IN  FILTER+2 


IOA  =  :?0404 

roe  =  :^o4os 

IOC  =  170400 

:ai  *  170402 

EALTi:  ”37  (R5)+ 

MOV  @<R5)+,@#ICC 
MOO  @<95>+,R0 

*av  POrRz 

oa°o:  mqi,.  .-*o  -  data'  re; 

ECS  R2*_CC30 
OOP::  MOV  RO,  RE 

CG°Z:  31”  7200 r  J#IGC 

C3  r  Q  I  ~  Q  “2  ** 

M0'-  0#IGIrR3 
MOO  R2 , 04 1 OA 
SUE  DA”A(RZ),R3 
MOV  .92, 9#  103 
ASR  R2 
ASR  R3 

ADD  R3,DA“A(RZ) 
DEC  R2 

SC 3  R2 , LC0P2 
:?R  LOG  PI 

aa?4:  rts  pc 

A”A  1  .  3LK**  1000 

.  END 


"ADD  OF  CH.  A  D/A 
JADD  OF  CH.  B  D/A 
‘  ADD  OF  CONTROL  3,  STA”US  REG. 

;add  OF  A/D  OUTPUT 

’ SETS  CH  NO  AND  EXT  TRIG 
;,?0  =  NO  OF  WE  I G-  ”5 

1 SE”  WEIGHT  VECTOR  ”C  ZERO 

JIN  IT  COUNTER 
' ”ES”  DONE  31” 

READ  A/D 

rZUTPIJ”  RAHP  ”C  CH.  A  D/a  “C  CHECa  SYNCH 
;R3  =  R3~«.EIGH7<R2> 

;”IL”ER  OUTPU”  TC  CH.  3  D/A 
J  SH IF”  R3  RIGhT  2  ”IM ES 

' ADJUST  WEIGHT 

’  DO  THIS  FOR  NUMBER  C-  weig-ts 

JGO  3AC>  and  S”AR”  A”  beginning  cr  vector 

",  RESERVE  RCCh  ”CF  WEIGHTS 
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STORE  PRESENT  POINTER  IN  DATA  MEM  5 
LP  SACL  3 

READ  IN  OLD  WEIGHT 
TPLR  4 

READ  A/D *  PUT  IN  DATA  MEM  a 
IN  a.  2 

THEN  INTO  ACC 
LAC  o*0 

CONVERT  FROM  OFFSET  BINARY  TO  2'S  COMPLEMENT 
XOR  7 

SUBTRACT  WEIGHT 
SUB  4.0 

STORE  RESULT  IN  DATA  MEM  6 
SACL  6 

CONVERT  BACK  TO  OFFSET  BINARY 
XOR  7 

sacl  e 

OUTPUT  TO  D/A 
OUT  8.2 

LOAD  HIGH  ACT  WITH  .2?  OF  OUTPUT 
LAO  a.  1.4 

NEW  WEIGHT  =  OLD  WEIGHT  f  .23  OUTPUT 
An  DP  * 

3ACH 

LOAD  ACC  WITH  LOCATION  OP  WEIGHT  IN  EXT  MEM 
L  A  C  0  •  0 

STORE  NEW  WEIGHT  IN  EXT  MEM 
TBLU  4 

INCREMENT  ACC  TO  POINT  TO  NEXT.  WEIGHT 
ADDS  1 

PUT  AUX  REG  0  OUT  TO  PORT  3  TO  CHECK  SYNCHRONIZATION 
SAR  0.9 
OUT  9*3 

CHCCN  IF  AT  END  OP  WEIGHT  VECTOR 
BANG  SINT 

I"  SO*  START  AT  BEGINNING 
L AR  0*3 
LAC  2*0 
INT  LA RP  1 

INSERT  DELAY  TO  PREVENT  DOUBLE  TRIGGER 
LARK  I  *  l SO 
EL  NOP 

BANZ  DEL 
L ARP'  0 

ENABLE  INTERRUPT 
El  NT 

GO  BACK  AND  WAIT  FOR  ANOTHER 
RET 
END 


Enclosure  2  to  Appendix  D 


*  NOTCH  FILTER  PROGRAM 

*  B.  b.  PETERSON-  IV  JUN  04 

.♦  TRIGGER  BY  EXTERNAL  INTERRUPT 

*  *  Of  UEIGHTS-1  MUST  BE  PUT  IN  DATA  MEM  3 

B  INIT 
NOP 

*  INTERRUPT  SERVICE 

B  ILP 

*  INITIALIZE 
INIT  LACK  0 

*  PUT  0  IN  DATA  MEM  0  TO  USF  ZEROING  TABLE  BELOW 

SACL  0 

*  SET  UP  A/H  AND  Li/A  CONTROL  REGISTER 

LACK  7 
SACL  10 
OUT  10 »w 

«  pur  IOOOOOOOO^opoOOO  CM  DATA  MEM  7  FOR  OFFSEF  BINARY  TO  2'S  COMPLEMENT  CON' 

! .  A  l !  h  ! 

•A  Cl.  ! 

i.Al.  l  •  i 0 

buCL  7 

*  mr  :  in  data  mem  i.  increment  in  table 

L  A  f.  K  r 
L-ACi  L 

»  PUT  LOCATION  OF  START  OF  TABLE  IN  DATA  MEM  1.'  •_  ■ 

LACK  251 
SACL  2 

*  SET  OVERFLOW  MODE 

ROOM 
LAPP  0 

*  ZERO  I A  I'LL 

*  Pin  STARTING  ADD  OF  TABLE  IN  ACCUMULATOR 

LAC  7  . 0 

*  PUT  I  OF  WEIGHTS- 1  IN  AUX  REG  0 

LAR  0.3 

*  POT  0  IN  THE  EXT  MEM  POINTED  TO  BY  THE  ACC. 

ZTB  TBLW  0 

*  INCREMENT  THE  ACC  BY  2 

ADDS  1 

*  DEC  AUX  REG  0  AND  LOOP  AGAIN  IF  NOT  ZERO 

BANZ  7TB 

*  MAIN  LOOT-  TO  START  AT  BEGINNING  OF  WEIGHT  VECTOR 

ML P  LAR  0 « 3 

LAC  2.0 

*  ENABLE  INTERRUPT 

El  NT 

*  WAIT  FOR  INTERRUPT 
WAIT  NOP 

B 


WAIT 
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ioao 

1085 

1060 

c 

c 

c 

1090 

1100 

1200 


1300 


XTSAflP 


loop i : 


conv: 


DO  1085  J  »  1 * NW 
Jl3 J  +  NW 
J2* J+2*NW 
J33 J+3*NW 
J4*  J-*-4*NW 
J5* J+5*NW 
J 6  3  J +6*N  W 

IF  (  I  FLAG  .  £Q  .  1  >  GO  TO  1080 

WRITE  (6. 1060)  J,  IX(  J  >  ,  IX(  Jl  >  ,  IX*  J2>  .  IX'.  J3)  *  IX<  J4)  ,  ix<  J5)  ,  IX(  J6> 
GO  TO  1085 

WRITE  (7,  1060)  J  *  IX  (  J  )  »  IX  (  J  1 ;  *  IX  i  J2 )  *  IX  (  J3  ;  *  IX  t  J4  >  ,  i;<  <  J5 )  *  IX  (  J6 ) 

CONTINUE 

FORMAT  (817) 

CALCULATE  AND  PRINT  OUT  SIGNAL  TO  NOISE  PATIO 
SNR»saiN/sauT 

WRITE<7, 1 100)  SGINrSOUT ,SNP 
FQPMA'  ■ 2F 1 0 . 2 ) 

WRITE  (7.1 200 ) 

FOPMAT  .  '  NEW  DATA'?  1  3  YES  -  0»n0  '  J 

READ  (5*140)  IFLAG2 

GO  TO  50 

STOP 

END 


Assembly  Language  Sampling  subroutine  called  from  main  proqram 


. GLOBL  XTSAMP 

B.  B.  PETERSON  *  17  NOV  83 

SUBROUTINE  TO  TAKE  A  GIVEN  NUMBER  OF  SAMPLES  FROM  THE  DTZ7S5 
ANALOG  I/O  BOARD  USING  EXTERNAL  TRIGGERING  AT  A  MAaIMUM  PATE  OF 
25*000  SAMP/S.  THE  SUBROUTINE  15  CALLED  BY  "CALu  XTSAMP ( N , NC *  I X ) ” 
WHERE  N*  #  OF  SAMPLES* 

NC  3  256*CHANNEL  #  >16.  AND  IX  IS  THE  DATA  VECTOR. 


IOC3 170400 
1 0 1  3 1 70402 
:  TST  iR5)+ 

MOV  0<R5)+*RO 
MOV  0 ( R5 )  +  *  0# IOC 
MOV  ( R5  )  •*■  *  R3 
BIT  #200*9# IOC 
BEG  LOOP  1 
MOV  0#IOI * R 1 
BIT  #200  * 0# I OC 
BEG  CONV 


' LOCAT I  ON  OF  CONTROL  REGISTER 
! LOCATION  OF  A/D  OUTPUT  REGISTER 

!R0«#  OF  SAMPLES 

ISET  CH#  AND  EXTERNAL  TRIGGER  ENABLE 
IR33STAFTING  address  of  data  vector 

*  TEST  DONE  B I T 
;loop  IF  NOT  SE' 

JREAD  to  CLR  DONE  BIT* throw  AwAY  IS'  SAMPLE 
’ TEST  DONE  BIT 
'  LOOP  IF  NOT  SET 


MOV  0# 10 1  * ( R3 ) ♦  JREAD  A/D  AND  PUT  DATA  IN  MEMORY 

DEC  RO 

BNE  CONV  JCHECK  TO  SEE  IF  N  SAMPLES  HAVE  BEEN  TAKEN 

RTS  PC 

.END 


in* . 
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250  CONTINUE 

C 

C  GET" N  SAMPLES  FROM  SELECTED  CHANNEL  AND  RETURN  AS  IX 

C 

300  CALL  XTSAMP(N,NC, IX) 

WRITE  (7,400) 

400  FORMAT  ('  SAMPLING  FINISHED') 

C 

C  INITIALIZE  PLOTTER 

C 

CALL  GON 
CALL  CLEAR 

CALL  PLIMITS  (.30,9.5, .3.6. 5) 

CALL  LI MI TS ( -FLOAT ( N ) / 8 . 0 , FLOAT ( N ) , - 1000 . 0 , 3000 . 0 ) 

500  NP* i N/NW ) 

CALL  MOVE  • -FLOAT(N) /9. 0,50.0) 

CALL  LABEL  ('  INPUT  ' , 2 , 0 , 0 ) 

C 

C  PLOT  INPUT 

C 

DO  525  1*1 ,N,NO 

CALL  LINE  ( FLOAT ( I )  , FLOAT 0 IX (  I  >  )  ) 

525  CONTINUE 

CALL  MOVE  (-FLOAT. N)/S. 0,2050.0) 

CALL  LABEL  ('  OUTPUT  2 , 0 , 0 ) 

C 

C  CALCULATE  AND  PLOT  OUTPUT 

C 

DO  800  1*1  ,NP 

DO  700  J=1 ,NW 

NS  *  (  I  —  1  )  -*NW+ J 
X*FLQAT i IX (NS) > 

Y  =  X-W<  J ) 

W<  J  >  *W<  J  )  -*-Y*G 
C 

C  PLOT  EVERY  NP  ' <rH .  POINT 

C 

IF  (ND1.LT.ND)  GO  T0  575 
CALL  LINE  ( FLOAT ( NS  > , Y+2000 . ) 

ND 1  *0 

575  ND 1 *ND 1+1 

C 

C  CALCULATE  INPUT  AND  OUTPUT  POWER  FOR  LAST  TWO  PERIODS 

C 

IF  (I.LT.NP-Z)  GO  TO  SOO 
SQIN*SQIN+ ■ 5*X*X/NW 
SOUT  =  SQUT  + . 5*Y*Y/NW 

700  CONTINUE 

800  CONTINUE 

CALL  GOFF 

IF  ( I FLAG . EQ . 0  >  GO  TO  1090- 
C 

C  PRINT  OUT  7  PERIODS  OF  INPUT  DATA  IN  7  COLUMNS 

C 


n  n  n  n  n  n  n 
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PROGRAM  PERDIC 
B.  8.  PETERSON 
:3  MAY  1984 

PROGRAM  TAKES  GIVEN  NUMBER  OF  SAMPLES  OF  AUDIO  WAVEFORM  AND 
CALCULATES  HOW  MUCH  OF  WAVEFORM  CAN  BE  CANCELLED  BY  PERIODIC 
WAVEFORM  AT  FUNDAMENTAL  FREQUENCY  OF  ENGINE 

arrays: 

IX  SIGNAL  INPUT 

C  COMMUNICATIONS  VECTOR  FOR  PLOTTING 

W  WEIGHT  VECTOR 

INTEGER  IX (10240? 

COMMON  C ( 20 ) 

DIMENSION  W( 312) 

SET  FLAG  TO  GET  DATA  FIRST  TIME 
IFLAG2  *  1 
ENTER  PARAMETERS 
WRITE  (7,100) 

FORMAT  ('  ENTER  #  SAMP,  #  SAMP/PLT,  #  WGHTS ,  CH  NO,  ADAP  GAIN') 
PEAD  (3,200)  N,ND,NW,nC,G 

SELECT  PRINTER  AND  PLOTTER  OPTIONS 

WRITE  (7,123) 

FORMAT  ('  0PTI0NS:0=N0  PRINTOUT , 1 =CRT , 2 3 PR INTER ' ? 

READ  (5,140)  I FLAG 
FORMAT  v 1 3 ) 

WRITE  (7,143) 

FORMAT  v'  0* PLOTS  ON  PLOTTER , 1 *CRT  '  ? 

READ  (5,140)  ID 
C( 10)=FL0AT( ID? 

ZERO  WEIGHT  VECTOR  AND  INPUT  AND  OUTPUT  POWER 

DO  130  J*  1  , NW 

W<  J ) =0.0 

SQIN»0 . 0 
SOUT -=0 . 0 

IF  (N.LT.O)  GO  TO  1300 
FORMAT  ( 4 1 G , F8 . 3 : 

SKIP  SAMPLING  TO  USE  OLD  DATA 

IF  ( I FLAG2 . EG . 0 >  GO  TO  300 
NC=236*NC+1G 

DELAY  BEFORE  SAMPLING 


DO  230  J* 1,30000 


UI;  XjiliJi'aj  j  |  I  j  int,i,n0. 


D8 


recorder.  The  plots  in  Figures  07  and  D8  were  produced  using 
this  program.  Enclosure  1  is  a  listing  of  the  program. 

The  filter  was  also  implemented  in  real  time  on  both  the  TMS 
320  and  the  LSI-11/2.  These  assembly  language  program  listings 
are  enclosures  2  and  3  respectively.  Because  the  TMS  320  has 
only  1 44  words  of  on-chip  data  memory,  it  was  necessary  to  access 
external  memory  using  the  Table  Read  and  Table  Write 
instructions.  These  instructions  require  the  address  to  be 
accessed  to  be  placed  in  the  low  order  accumulator.  This  results 
in  relatively  complicated  code  to  implement  a  very  simple 
algorithm.  Despite  this  inefficient  use  of  the  TMS  320,  the 
algorithm  can  still  run  at  sampling  rates  of  up  to  150-  kHz. 

The  MACRO-11  (PDP-11  assembly  language)  program  in  enclosure  3 
illustrates  the  simplicity  of  the  algorithm.  On  an  LSI-11/2  this 
program  will  run  at  sampling  rates  of  up  to  11  kHz. 

Cadet  1/e  FAVERO  is  presently  attempting  to  implement  the 
filter  on  an  AIM-65  microcomputer  at  sampling  rates  of  at  least  8 


kHz. 


30^  LOCUS  OF  S~H  ORSER  F 

REFERENCE  1NPUT  IS  R  SQURRE 

Figure  Da 


.  <  • 
✓n  n  v 


samples,  the  transfer  function  from  xl  (k)  to  y(k)  is  given  by: 


Y(z) 

z”  - 

1 

X,Cz) 

N 

z  + 

M-T 

Z 

a  dj  zj  -1 

where 

dJ 

SIaJ 

M 

j=o 

x2 

(i)  x2(i+j) 

i  =  0 

s  N  x  the  autocorr elat ion  of  x 

Figure  D5  is  a  root  locus  plot  as  the  adaptive  gain,  a,  is  varied 
for  N  =  3  and  for  x ^ C tc )  a  square  wave  and  Figure  D6  is  the 
frequency  response  of  the  same  filter.  Mow  the  width  of  each 
notch  is  a  function  the  amplitude  of  the  corresponding  harmonic 
in  x2(k)  in  addition  to  the  adaptive  gain.  This  suggests  the 
possibility  of  using  apriori  knowledge  of  which  harmonics  are 
strongest  in  the  noise  component  of  x^k)  to  specify  an  x2(k)  for 
more  effective  filtering. 

Another  possible  interesting  application  of  this  type  of 
filter  is  for  very  narrow  notch  and  bandpass  recursive  digital 
filters.  The  major  advantage  is  that  even  with  low  precision 
fixed  point  arithmetic  (8  bit,  for  example)  and  for  poles 
virtually  on  the  unit  circle,  the  filter  is  completely  stable. 

The  filter  has  been  implemented  in  three  different  ways. 

First,  u'Sing  off-line  processing  in  FORTRAN,  a  data  vector  is 
obtained  by  externally  triggering  the  LSI-11/2  A/D  converter  and 
then  the  filter  implemented  and  the  data  plotted  on  a  strip  chart 
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INITIALIZE  COEFFICIENTS  TO  CALCULATE  WEIGHT  ■,'EC"rQR 


AND  ZERO  WEIGHT 

VECTOR 

DO  I--1,NW2 

°2  (  T  '  r  « 

OE  - 10 

CR  (  I )  •»  1  . 

AC-i  ,*> 

C I  ( I  )  «  1 . 

OC-  1  13) 

XH ( I ) =0 . 

rv 

Vu ( I ? -o . 

CONTINUE 

N8L0CS*N/NW 

DO  1*0  J  -  1  r  N8L0CS 


JNW=  <  J  -  1  )  •"•MW 

FILL  '.'ECTPS  TQ  CALCULATE  CCT 

° I *0 . 0 
°C~0 . 0 

DO  120  I  - 1  ,  NW2 

I A 1 “  2*I*JNW 


I A  =  I A 1  - 1 

XI  (  I  )  - 

Cl  H(STf  i  i  |t(il  Z'lftSI 

Y 1  (  I  )  ^ 

ci  n^T/H  ('T'A1  >  n 

X2  (  I )  - 

ei_0AT ' 1 2 ( I A  ? / 205 ) 

Y2  <  I  ?  - 

Cl  0AT<  12  (  IA1  )  / 20? ) 

CONTINUE 

CALL  RFFT  (  NM  ,  NW  •  NNAX  »  X  1  ,  Y 1  ,  WR  ,  U I  1 . 0  ) 
CALL  RFCT ( NW  t NW , NNAX • X2 , Y2  ,  WR , Wl ,1.0) 
DO  130  I  -  1 , NW2 


CALCULATE  nc  QyTO'JT 


VQ  f  r  \  -  VI  (  r  ^  _  vu  f  r  t  v  2  r  I  I  '  ~  YZ  (  I  ' 

YQ  (  l  \  s  y  1  '  T  '  -  Yu  ^  I  )  *v2  f  I  '  —  VW  (  I  )  *  Y2  (  I 


CHECK  CCR  S  M  ALL  VAL 

UES  AND 

SE 

t  rn  ; 

t  rz ; 

~0 

PREVENT  UNDERFLOWS 

T  C 

(Aes <  XI  <  I  )  > . LT. 1 

#AC-1  |)  N 

v 1 

(  T)  -O 

# 

T  C 

'■  ASS  (  Y 1  (  :  )  )  .LT.  : 

4 AC- 1 A  \ 

Y 1 

/  ;  >  «?0 

4  (  _*> 

I F 

<  ASS ( Y2 ( I )  )  lt . i 

.  AE- 1 0 > 

VC 

(  I  >  *=0 

#  A 

I F 

( ASS < X2 ( I ) > .Lt. 1 

•  0  E  ”  1 0 ) 

v  r* 

r  r  >  -  0 

#  0 

CALCULATE  INP'JT  and 

OUTOL,r’ 

pp 

'•IS3 

SI:DI+Y1  (  T  )*v’  ,I'-*-YlfI'J‘vlfI' 

pg  =  pg  +  vg (  r  ;*VQ(  T  \  * YO  (  I  )  *Y0  I  ) 


update  weights 


P2  i  I  )  =  Qf*,G-K-D2'’  I  )  *  G-*’  (  X  2  '  I  '  •»  v2 '  I  )  Y2  '  I  )  *  v 
C.9  <  I  )  s0NG*Cr?  {  I )  ■+■□■*  ( y  1  <  I  )  *X2  '  I  )  -*■ y  1  (  I  )  •»  Y 
Cl  (  I  >  =ONG*C I  <  I  ?  +G-*  < M2  (  I  )  *Y1  (  I  )-Xl  <  I  )*Y 
XH ( I ) =CR ( I ) / °2 (  I  ) 


tin  ri 
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FFT  subroutines  called  frcm  main  program 


cct  SUSPOUT INES 

o  o  ,  qctcdcq^i  =■  era 

subroutine  ?cc"r  •: m  .  m  , nc*ax  ,  x i  y /  nr f  w i ,  d:  > 

SUBROUTINE  TO  CACLULATE  N  POINT  PEAL  USING 

N/2  °C  I  NT  COMPLEX  F-T 

REF  I  "TMS  020  DIGITAL  SIGNAL  PROCESSING  SEN I  NAP  NOTES. 
TEXAS  INSTRUMENTS,  INC.,  1SB3. 

N*  #  OF  PEAL  DATA  POINTS  =  2**M 

X  IS  REAL  PART  qf  F^t  <n/2  VECTOR )  Y  IS  IMAG  PA PT 
WR  AND  MI  APE  VECTORS  OF  COSINES  AND  SINES  OF  =122  NMA 
D I =  1.  FOR  FORWARD  FFT ,  -i.  FOR  INVERSE. 

FOR  FORWARD  FFT  X  AND  Y  SHOULD  BE  GENERATED  BY 
X ( I  >  =  DATA <2*1-1  )  AND  Yd)  =  DAT 6(2+1) 

WHERE  DATA( )  IS  THE  INPUT  VECTOR  OF  REAL  POINTS 

REAL  X< 1 ? , Y< 1 > ,WR( 1 ) ,WI ( 1  ) 

N2=N/2 
NM1  -  N— 1 

IF  (DI.LT.O.O)  GO  t q  in 

CALL  ==-T2<N2,HNl  ,NNAX,X.,  Y,WR,UI  ,  1 .0). 

NMON*NMAX/N 

N4=N2/0 

Y ( N2  + 1 ) =Y( 1 ) 

X  <  N2  +  1 ) =X ( 1 ) 

I  Os  i 

DO  20  I  =  1  ,  'Ji+S 

NHI -N2-I*2 
9a= v '  i  >  *y. Nr*  I  ) 

XA  =  Y<  I  -Y<NMI 
R9  -  (  X  ( NM  I  >  —  x  <  i  >  \ 

'/o  -  (  _y  (  I  >  -  v  (  NM  I  ) 

£pusUJft  (TO) 

SPU=-W I (IP) 

IF  (DI.GT.O.O?  GO  T0  If 
RS=-«9 

xs»-xs 

5pu-_ spu 
IPs  I P+NMON 
RC  -SPI-IJ*-RB-CPUJ‘X3 
XC  -SPH*XB*CPU'»RS 
X(I)*.5*(RA+RC? 

Yd.)  =  .f~(XA-»-XC) 

X ( NM I )  =  . f  * ( RA-RC ) 

Y ( NM I ) = . f * ( XC-XA ) 

CONTINUE 

IF  (DI.GT.O.O)  GO  TO  1000 

CALL  FFT2 ( N2 , MM  1 , NMAX, X , Y , HR, W! .-1.0) 

RETURN 


nnnnnnnnnnnnn 
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FFT  subroutines  called  frcn  main  program 


COMPLEX  F^T  WITH  TABLE  LOOKUP 

PEFIC.  S.  BURRIS  AND  W.  PARKS  /  ”DFT/FffT  AND  CONVOLUTION 
ALGORITHMS , "  JOHN  WILEY  L  SONS,  NEW  YORK,  ISO* 

CALLING  VARIABLES: 

N  NUMBER  OF  COMPLEX  POINTS 

M  L0G2  OF  N 

NMAX  MAXIMUM  VALUE  OF  N  FOR  INTERPRETING  LOOKUP  TABLE 
X  L  Y  REAL  AND  IMAGINARY  PARTS  OF  FFT 
WR  TABLE  OF  COSINES 

WI  T ABLES  OF  -SINES 

DI  *1.0  FOR  FORWARD  PFT,  -1.0  FOR  INVERSE 


40 

50 

100 


101 

102 


SUBROUTINE  FFT2  <  N , M , NMAX , X , Y , WR , W I , DI  ) 
REAL  Y  <  1  )  ,  X  <  1  )  ,  WR  <  1  1  ,  W I  >.  1  > 

NMON*NMAX/N 
N2 -N 

DO  ICO  K  *  1  »  M 
Ml  -N2 

M->  _  —  /  — 

I E  *  M  /  N  1 
IA=  1 

DO  50  J»1,N2 

I A  1  -  (  I  A—  1  >  •*NMQN*i 
*  C-MR(IAl) 

*  S*-DI*WI CIA1) 

DO  40  I-J,N,N1 
L - I*N2 

«r- y  (  r  \  _v  ( t_  \ 

y ( r I ) *  X ( L > 
yr*y  (  r  \  _y  ( t_  \ 

Y  (  I  ?  a  Y  C  I  >  *  Y  (  L  ) 

Y  (  L  )  =  C*Y"r-S»v'"r 

CONTINUE 

CONTINUE 

CONTINUE 

J*1 

N 1 *N— 1 

DO  10^  1*1 , N 1 

IF  ( I . GE. J  >  GO  T0  101 
XT“Y ( J ) 

X( J ) *X<  I  ) 

Y  (  T  \  -  y  T 
VT - Y ( J ) 

Y  '  J  '  a  Y  <  I  ) 

Y  (  t  ^  - v-r 

K -N/2 

!=■  'K.GE.J?  GO  rG  100 
j  -  J-K 

y  9  W  /  *» 

GO  T0  102 

J  *  J  ♦  ¥ 


103 
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CONTINUE 

IF  (DI.GT.O.O)  GO  tq  io6 
DO  105  I  -  1  .•  Nl 

y.  <  i )  *  x  ( i )  /  m 

Y(I>*Y(I)/N 

CONTINUE 

RETURN 

END 

SUBROUTINE  PACK  <  N ,  NMAX  .  IX  ,  X ,  Y  ,  IN  ,  WIN > 

PUTS  INTEGER  DATA  IN  REAL  ARRAYS  SU I TABLE  FOR  RFc"r 
MULTIPLIES  EACH  DATA  POINT  BY  HAMMING  WINDOW 

REAL  X  (  1  )  ,Y(  1 >  ,  WIN*  1 ) 

INTEGER  IXC) 

NMON*NMAX/N 

N2“N/2 

DO  10  I = 1 , N2 

-"2-1  *2 

x7l)=PL0AT  <  IX  (12-in 
Y  (  I  ) =eLQAT  (IX (12)) 

IF  (IW.EQ.O)  GO  T0  10 
1 2  =  <  12-2)  *NM0N+1 
X( I )=X( I )*WIM( 12) 

Y ( I ) * Y (  I ) *WIN<  MMOM+ 13 ) 

CONTINUE 

RETURN 


SU8  ABOVE 


Appendix  F  LI5S  Adaptive  Filter  FORTRAN  Program  v,±iere  Audio  Signal  is  added 


C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

50 

100 

200 


200 

325 


350 


370 

450 

500 

S00 

C 

C 


in  Software 

P90GRAM  ADAPT 
B.  8.'  PETERSON 

PROGRAM  READS  IN  ANALOG  DATA  FROM  AN  F*  TUNER 

DECK  USING  SUBROUTINE  XTSMP2.  IT  THEN  SAMPLES  TWO  CHANNELS  F»OM 
A  TAPE  DECK.  THE  FM  TUNER  SIGNAL  IS  MULT.  BY  A  GAIN  AND  ADDED 
TO  ONE  NOISE  CHANNEL  TO  GENERATE  THE  PRIMARY  INPUT. THE  OTHER  NOISE 
CHANNEL  IS  ADAPTIVELY  FILTERED  AND  SUBTRACTED  FROM  THE  DELAYED 
PRIMARY  INPUT.  THE  WEIGHTS  IN  THE  ADAPTIVE  FILTER  ARE  ADJUSTED  SO 
AS  TO  MINIMIZE  TH£  OUTPUT  POWER.  THE  REFERENCE  INPUT,  THE  PRIMARY 
INPUT,  THE  FM  TUNER  SIGNAL  AND  THE  OUTPUT  ARE  JLOTTED  ON  THE  CRT, 
THE  HP  7470  PLOTTER  OR  BOTH.  THE  SIGNAL  TO  NOISE  RATIOS  OF  BOTH 
INPUT  AND  OUTPUT  ARE  CALCULTED  AFTER  Twg  WEIGHTS  HAVE  HAD  SOME 
TIME  TO  CONVERGE. 

arrays: 

IY  PRIMARY  INPUT 

IX  REFERENCE  INPUT 

ISIG  INTEGER  SIGNAL  VECTOR 

H  WEIGHT  VECTOR 

c  communications  vector  for  °ltl:b 

INTEGER  IX < 4000 > , IY< 4000) , ISIG <4000) 

REAL  H  < 100) 

COMMON  C ( 20 ) 

INPUT  PARAMETERS  * 

WRITE  < 7 , 100 ) 

FORMAT  ('  NEW  DATA?,  RESET  WEIGHT  VECTOR?  i*YES,0«NO'> 

READ  (5,200)  IF1 , IF2 
FORMAT  (2IS ) 

IF  (IF1.EQ.0  .AND.  IF2.EQ.0)  GO  TO  450 
WRITE  <7,30 0? 

FORMAT  ('  #  SAMP, REF  !NP,PRI  INP,SIG  INP-#  WTS  '  ) 

READ  (5,325)  N , NC 1 , NC2 , NC3 , NW 
FORMAT  <5IG> 

NT  aNW/2 
N2=N/2 

WRITE  (7,350) 

FORMAT  ('  PLOT  RAW  DATA*?  ,  0  =  ?4?0 , 1  =G  I  Ci  I  ,  2  =  BOTH  ,  -  :  =M0  PLOT, STEP') 
READ  (5,200)  IC10,ISTP 
STP=FLOAT( ISTP) 

IF  (IC10.LT.0)  GO  tq  370 
C(10)*FL0AT  (IC10) 

NC!*25S*NC1+1S 
NC2=25S*NC2+ IS 
NC3=256*NC3+1S 
WRITE  (7,500) 

FORMAT  ('  ADAPTIVE  GAIN  <E  cORMAT) , SIGNAL  GAIN  < I  FORMAT)') 

READ  (5, BOO)  G,IG 
FORMAT  (El  4. 4, IS) 

IF  (IF2.EG.0)  GO  TO  300 

ZERO  WEIGHT  VECTOR 


i  ibi  i  Tb 


FI 


F2 


C 

DO  700  I=1,NW 
700  H ( I >  =0 . 0 

900  IF  (IF1.EQ.0)  GO  TC  325 

C 

C  DELAY  8EF0RE  SAMPLING 

C 

DO  850  1=1,10000 
850  CONTINUE 

C 

C  SAMPLE  SIGNAL  FROM  TUNER 

C 

CALL  XTSMP2<N,NC2,NC2, ISIG, IY) 

C 

C  SAMPLE  TUO  CHANNELS  FROM  TAPE  DECK 

C 

CALL  XTSMP2(N,NCl ,NC2, IX, IY) 

825  IF  (IC10.LT.0)  GO  TO  SOO 

CALL  PLIMITS  (0.0,10.0,0.0,7.0) 

CALL  LIMITS  ' 0 . 0 , FLOAT < N ) , - . 3*STP , 3 . S*STP ) 

C  IF  (IF1.E0.0)  GO  TO  900 

DO  875  1  =  1  ,N 

CALL  LINE  ( FLOAT ( I ) , FLOAT ( IX (  I  )  )  ) 

875  CONTINUE 

CALL  MOVE  ( 1 . 0 , S  T  P  J 

DO  890  I  *  1  , N 

ISIG(I>  =  IG*ISIGd> 
i  y  ( i- )  =  i  y  ( :  '+is:g(  i  > 

CALL  LINE  '■  FLO  A”  (  I  )  ,  FLOAT  (  I STD  + 1 Y  ( I )  )  > 

390  CONTINUE 

CALL  MOVE  ■ 1.0,2*570? 

QFCSE” =  2.0*ST3 
DO  895  I  =  1  ,  N 

CALL  LINE  !  CLQA~  (  I  )  .•  FLOAT  < IS IG (  I )  >  +G- rSE~  ) 

395  CCN'INUE 

pprcr--T  _  n»c-3 

CAL-  MOVE  (  :  .0 • Ce-=E”> 

C 

2  22^C  ■.TLA--  3L  -  NO :  SE ,  OUTPUT  NOISE,  AND  DESIRED  OUTPUT. 

C 

900  SGIN=0.0 

SQOUT  =  0.0 
SOS IG=0 . 0 

c 

C  MAIN  ADAPTIVE  FILTER  LOOP 

C 

DO  1000  I*NUI+1,N 
FOUT  =  0.0 
C 

C  CALCULATE  TL^ER  OUTPUT 

C 

DO  950  J=1 ,NW 

IS=< I-J) 

950  FOUT  =  pOUT +H ( J)*IX< IS) 


ID® ( I— NT ) 

OUT®  I Y ( ID ) -f OUT 


C 

C 

C 

C 


C 

C 

c 

980 

C 

C 

c 

980 

975 

1000 


1050- 

1099 

1100 


SQUARE  AND  LON  PASS  PILFER  INPUT  AND  GUTPUT 
FOR  LAST  HALF  OF  DATA 

IF  ( I .LT.N2)  GO  TO  geo 

SGQUT=.98  +  SGQUT+.02-*<QUT-ffLQAT<  ISIG< ID) >  ) 

R IX  =  FLOAT  < IY( ID)-ISIG( ID) > 
SQIN=,98*SGIN-»-.02*RIM*RIX 

SGSIG®.98*SQSIG+.02*FL0AT( I S I G  ( I ) >*FLQAT( ISIG 

SNRIN=SGSIG/SQIN 

SNROUT  =  SQSIG/SQOUT 

P'-0T  OUTPUT 

CALL  LINE  < FLOAT (  I  > r  OUT +0FC5ET > 

AD APT  WEIGHT  '.'ECTOR 

DO  975  J®  1 , MW 

IS®< I-J  » 

H ( J ) ®u<  j )>g*OUt*IX( IS) 

CONTINUE 
CALL  TURNOF 

WRITE  (7.1050)  SNRIN, SNRQL'T 
FORMAT  (2F10.5) 

GO  TO  50 

STOP 

END 
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